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Abstract 

The  remote  observations  of  the  temporal  and  spectral  characteristics  of  the 
infrared  (IR)  emissions  from  exploding  ordnance  have  been  correlated  with  explosion 
conditions.  A  Bomem  MR- 154  Fourier  Transform  Interferometer  with  two  detectors, 
InSb  and  HgCdTe,  was  used  to  record  spectra  in  the  1.3-20  pm  range.  Data  was 
collected  at  spectral  resolutions  of  16  cm'1  and  4  cm'1,  and  temporal  resolutions  of  0.045 
s  and  0.123  s,  respectively.  The  data  files  range  in  size  from  900  Kilobytes  to  several 
Megabytes.  These  are  reduced  to  2-dimensional  representations  of  temporal  features  that 
are  less  than  100  Kilobytes. 

The  data  collected  from  efforts  indicate  the  possibility  of  characterizing  event 
species  through  one  or  more  derived  temporal  features.  Each  event  data  matrix  contains 
three  dimensions  of  information  describing  radiance  as  a  function  of  frequency  and  time. 
The  observed  data  is  first  corrected  for  atmospheric  losses  to  convert  apparent  radiance  to 
emitted  radiance.  The  data  is  then  adjusted  to  remove  background  radiance.  Finally,  the 
corrected  data  is  fit  to  a  Planckian  distribution  function  to  compute  event  temperature  and 
fractional  field  of  view.  Temporal  profiles  of  temperature  and  fractional  field  of  view  are 
created  to  describe  each  event.  The  temporal  profiles  of  each  explosive  type  are 
compared  to  other  explosive  types.  Certain  explosive  types  indicate  an  afterbum  feature 
on  their  temperature  profiles.  The  afterburn  feature  wasn’t  apparent  on  the  temperature 
profiles  of  other  types.  Additionally,  the  temporal  evolution  of  fractional  field  of  view 
was  unique  for  each  explosive  type. 


I.  Introduction 


Overview 

NASA  defines  remote  sensing  as,  “The  acquisition  and  measurement  of 
data/information  on  some  property(ies)  of  a  phenomenon,  object,  or  material  by  a 
recording  device  not  in  physical,  intimate  contact  with  the  feature(s)  under  surveillance; 
techniques  involve  amassing  knowledge  pertinent  to  environments  by  measuring  force 
fields,  electromagnetic  radiation,  or  acoustic  energy  employing  cameras,  radiometers  and 
scanners,  lasers,  radio  frequency  receivers,  radar  systems,  sonar,  thermal  devices, 
seismographs,  magnetometers,  gravimeters,  scintillometers,  and  other  instruments.”  (22, 
2-4).  By  this  definition,  remote  sensing  can  be  as  simple  as  the  human  eye  sensing  the 
reflected  light  from  this  text  as  it  is  being  read  or  as  complicated  as  satellite 
measurements  of  distant  parts  of  the  universe.  As  technology  has  grown,  so  have  the 
opportunities  for  ways  to  use  remote  sensing.  No  longer  is  remote  sensing  limited  to  the 
visual  portion  of  the  electromagnetic  spectrum.  Today  radiation  can  be  captured  and 
recorded  from  short  wavelength  x-rays,  through  visual  and  infrared,  and  beyond  to  longer 
wavelength  radio  waves.  Each  spectral  range  can  be  used  to  detect  unique  information 
about  unknown  objects  or  events.  For  example:  X-rays  detect  solar  flares,  visual 
photographs  provide  military  intelligence  or  pictures  for  disaster  relief,  infrared 
determines  agricultural  crop  conditions  as  well  as  theater  missile  defense,  and  radio 
waves  monitor  changing  conditions  in  the  polar  icecaps.  Remote  sensing  has  been  used 
in  various  applications  including  medical  diagnostics,  weather  satellite  imagery,  terrain 
mapping  of  Mars,  to  identifying  groceries  in  the  supermarket  using  barcode  scanners. 
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Remote  sensing  is  widely  used  in  the  Department  of  Defense  (DOD).  One  of  the 
first  uses  of  remote  sensing  by  the  US  Army  Air  Corps  was  photographing  troop 
movements  in  WWI.  The  DOD  now  uses  sensing  techniques  for  battlefield  management, 
battle-space  characterization,  weapons  guidance,  technical  intelligence  and  threat 
identification.  Today,  remote  sensing  is  not  limited  to  visual  pictures  of  troop 
movements,  but  includes  sensing  most  of  the  electromagnetic  spectrum.  Certain  regions 
of  the  spectrum  provide  information  about  distinct  events  not  apparent  in  other  regions  of 
the  spectrum. 

This  thesis  continues  the  work  of  Orson  (17),  “Collection  of  Detonation 
signatures  and  Characterization  of  Spectral  Features”.  His  research  focused  on  the 
infrared  collection  of  electromagnetic  signatures  from  fifty-six  detonation  events.  These 
events  were  collected  over  the  spectral  range  of  500  to  6000  wavenumbers  (cm1)  with 
.047  seconds  time  resolution.  Collection  was  accomplished  with  a  Bomem  MR-154 
Fourier  Transform  Infrared  spectroradiometer.  These  events  were  a  combination  of 
dynamic  F-18  delivered  ordnance  as  well  as  static  detonations.  The  research  was  part  of 
infrared  Measurement  and  Signals  Intelligence  (MASINT)  detection  of  explosions.  This 
MASNIT  research  is  part  of  a  cooperative  effort  between  the  Navy  Tactical  Exploitation 
of  National  Capabilities  (TENCAP)  office  and  the  National  Air  Intelligence  Center’s 
Cobra  Brass  program  office  (NAIC/DXDI).  This  program  has  sponsored  a  series  of  tests 
where  ordinance  was  detonated  under  field  conditions.  These  efforts  explore  the 
potential  use  of  space-based  MASINT  sensors  to  support  combat  units  in  the  area  of 
Battle-Space  Characterization. 
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The  source  of  detonation  data  is  from  munitions  detonations  performed  during  the 
Navy  TENCAP  sponsored  testing  of  munitions.  This  series  of  tests  was  named  Radiant 
Brass.  This  series  consisted  of  five  tests  conducted  at  the  Navel  Air  Station,  (NAS) 
Fallon,  NV.  Testing  started  in  June  1998  and  was  completed  October  1999.  Throughout 
the  program  various  contractors  and  government  agencies  have  supported  Radiant  Brass 
by  deploying  sensors  to  the  test  sight  collecting  event  signatures.  The  primary  data 
requirements  for  these  collection  teams  were  radiometric  in  nature.  For  these  tests, 
spectral  information  was  a  secondary  objective. 

Blackbody  Radiation  and  Atmospheric  Transmittance 

A  unique  opportunity  for  remote  sensing  exists  in  the  InfraRed  (IR)  portion  of  the 
electromagnetic  spectrum.  Within  the  IR  portion  of  the  spectrum  are  the  peak  intensities 
of  radiation  from  emitters  including  300  Kelvin  to  thousands  of  degrees  Kelvin.  The 
temperature  of  a  graybody  can  be  related  to  the  spectral  radiation  emitted  by  an  event  via 
Equation  1.1,  the  Planck  Distribution: 

27Thc2  1 

S(X)  =  s(X)^  1  (1.1) 

A  e  /m  _  \ 

where 

S(A)  =  Spectral  Radiance  as  a  function  of  wavelength 

e(X)  =  emissivity  as  a  function  of  wavelength 

X  =  radiation  wavelength 

h  =  Planck’s  constant  =  6.626  *  10'34  W  sec2 

T  =  absolute  temperature  in  Kelvin 

c  =  velocity  of  light  =  2.9979  *  108  m/sec 

k  =  Boltzmann’s  constant  =1.38  *  10‘23  W  sec/K 

The  wavelength  of  maximum  emittance  (Xm)  in  microns  is  given  by, 

Xm  =  a/T  (1.2) 
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A  simplified  form  for  total  emittance  (M),  of  the  Planck  distribution  are  given  by  the 
Stefan-Boltzmann  Law: 

M  =  aT4  (1.3) 

where 

a  =  constant  =  5.669  *  1(T8  W  m2 K'4 


If  s(A)  =  1,  Equation  1.1  describes  the  radiation  emitted  from  a  black  body.  When 
£?±\,  but  is  constant,  the  output  from  Equation  1.1  is  a  called  a  ‘gray’  body.  When  s  *  1 
and  is  wavelength  dependant,  it  is  described  as  a  selective  radiator. 

Many  texts  and  examples  of  spectrum  use  wavelength  as  its  unit  of  length 
measurement.  In  spectroscopic  discussions,  frequency  in  the  form  of  wavenumber  is 
more  commonly  used.  The  conversion  between  wavenumber,  v,  and  wavelength,  A,  is: 


X(fim ) 


10000  ^ 

v(cm~l)  j 


(1.4) 


Since  the  spectroradiometer  software  outputs  all  information  in  cm'1,  these  units 
will  be  used  primarily  in  this  thesis.  It  is  important  to  note  that  the  wavelength  to 
wavenumber  conversion  is  not  linear  in  wavelength.  This  changes  the  shape  of  the 
Planckian  curve.  Higher  wavenumber  equates  to  lower  wavelength  and  the  wavelength 
unit  of  measure  trends  in  the  opposite  direction.  The  Planck  Distribution  function  as  a 


function  of  wavenumber  becomes: 


S(V)  =  *(V) 


2hc2v 3 

chv/ 

e  /kT-\ 


(1.5) 
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Figure  1.1  shows  the  shapes  and  relative  intensity  of  spectral  radiance  of  black 
bodies  at  different  temperatures.  The  Figure  shows  four  temperatures;  500,  700,  900  and 


1 100  degrees  Kelvin. 
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Figure  1.1  Spectral  Radiance  vs.  wavenumber 


Also  within  the  IR  region  of  the  spectrum  there  are  “transparent”  areas  where  the 


atmosphere  does  not  significantly  absorb  electromagnetic  radiation  and  “opaque”  regions 


where  energy  at  a  given  frequency  will  be  largely  absorbed.  Phillips  Laboratory  Expert 


User  Software  (PLEXUS)  from  the  Mission  Research  Corporation  (MRC)  was  used  to 


create  atmospheric  transmittance  profiles.  PLEXUS  is  a  knowledge-based,  highly  user- 


oriented  software  platform  that  integrates  and  widens  the  accessibility  of  the  Air  Force 


Research  Laboratoiy  (AFRL)  family  of  atmospheric  and  astronomical  background  codes. 


PLEXUS  provides  an  intelligent,  intuitive,  graphical  environment  for  setting  up  and 


running  the  AFRL  Codes  and  analyzing  their  output.  The  AFRL  atmospheric  codes  are 


the  DoD  standard  models  for  computing  the  spectral  radiance  and  transmittance  in 


various  regions  of  the  atmosphere.  The  atmospheric  and  celestial  codes  in  PLEXUS  are 
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official  releases  by  AFRL.  MODerate  TRANsmission  (MODTRAN)  was  the  AFRL 
atmospheric  model  used  in  PLEXUS  for  this  research.  MODTRAN  code  models  the 
transmission  and  radiance  generally  applicable  below  an  altitude  of  50  km. 

Transmittance  is  used  to  describe  the  portion  of  source  electromagnetic  radiation  received 
at  the  detector  after  traveling  through  a  given  path.  Transmittance  is  a  unitless  number 
ranging  from  zero  to  one.  A  value  of  one  represents  no  signal  loss  and  a  value  of  zero 
represents  a  total  loss  of  signal  between  the  source  and  the  detector. 

Figure  1 .2  shows  a  PLEXUS/MODTRAN  transmittance  profile  of  the  atmosphere 
as  a  function  of  wavenumbers.  The  transmittance  profile  shown  in  Figure  1.2  represents 
the  portion  of  electromagnetic  energy  that  passes  between  a  sensor  and  a  detonation  event 
two  miles  away.  The  profile  shown  in  Figure  1 .2  was  created  using  the  atmospheric 
conditions  recorded  at  Fallon  NAS  during  testing.  The  profile  represents  a  4200  meter 
line  of  sight  path  through  the  atmosphere  10  meters  above  the  ground.  Figure  1.2  has  a 
spectral  resolution  of  4  cm'1.  The  frequencies  with  values  of  zero  are  where  the 
atmosphere  has  completely  absorbed  the  radiation  at  that  frequency.  The  primary 
absorbers  in  the  IR  are  H20  and  C02.  Figure  1.2  also  shows  which  molecular  species  is 
primarily  responsible  for  absorption  in  opaque  frequency  intervals.  There  are  also  areas 
where  the  atmospheric  transmittance  is  close  to  one.  In  these  frequencies,  the  radiation 
received  by  the  detector  is  in  the  atmospheric  and  spatial  conditions  used  is  nearly  the 
same  as  the  radiation  detected  if  there  were  no  atmosphere.  Theoretically,  these 
frequencies  represent  where  the  data  is  most  accurate. 
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Figure  1.2  Atmospheric  Transmittance  Profile 

By  multiplying  the  transmittance  profile  by  a  Planckian  blackbody  curve  at  a 
given  temperature,  an  expected  spectral  signal  can  be  created.  Consider  the  1100  Kelvin 
curve  from  Figure  1 . 1  traveling  through  the  atmosphere  to  our  detector.  Figure  1.3 
represents  what  an  1100  Kelvin  blackbody  would  be  for  the  atmospheric  and  spatial 
conditions  used  to  create  the  Figure  1.2  transmittance  profile.  The  detector  would 
intercept  this  signal  if  the  blackbody  completely  fills  the  detector’s  field  of  view. 
However,  the  detector  can  be  calibrated  and  amplify  the  signal  to  approximate  the 
expected  event  as  discussed  in  chapter  II.  Figure  1 .3  is  created  by  multiplying  the 
unit  less  numerical  value  of  the  transmittance  by  the  intensity  value  of  a  blackbody  at  a 
given  temperature  for  each  frequency  value. 
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Figure  1.3  1100  Kelvin  Intensity  with  Atmospheric  Transmittance 

This  predicted  signal  can  compared  to  actual  data  from  one  time-step  just  after  bomb 
tion  in  Figure  1 .4  shown  below. 
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Figure  1.4  Raw  Data  for  One  Time  Step 


Figure  1 .4  represents  the  spectral  signal  at  one  time  interval  within  our  total  data 
set.  In  this  case,  the  spectral  profile  at  a  time  interval  just  milliseconds  after  detonation  is 
shown  to  approximate  an  1 100  Kelvin  blackbody  profile.  The  goal  of  this  research  will  be 
to  interrogate  each  time  step  and  identify  a  pattern  for  the  entire  data  set.  To  do  this  it  is 
important  to  understand  what  is  expected  from  the  observed  event.  Since  this  thesis 
evaluates  the  bomb  detonation  data  from  the  Radiant  Brass  tests,  a  quick  discussion  on 
explosion  mechanistics  is  necessary. 

Explosion  Mechanistics 

A  detonation  produces  three  primary  ‘sensed’  features,  an  explosion  initiation,  a 
secondary  afterbum,  and  a  decay  period  (Cooper,  133).  The  explosion  initiation  is  caused 
by  the  nearly  instantaneous  conversion  of  the  explosive  from  chemical  potential  energy  to 
kinetic  energy.  The  time  scale  of  this  feature  is  less  than  1  ps. 

The  afterbum  feature  is  the  result  of  the  residual  products  of  the  detonation.  In 
an  under  oxidized  explosion  these  residuals  are  themselves  fuels.  After  the  initial 
pressure  pulse  of  the  explosion  reaches  equilibrium  pressure  with  the  surrounding 
environment,  oxygen  and  air  remix  with  the  residual  explosion  by-products.  This 
mixture  of  oxygen  and  residual  fuel  ignites  and  causes  the  energetic  afterburn  fireball. 

This  feature  starts  milliseconds  after  the  initiation  feature  and  can  last  for  up  to  a  second. 
The  third  feature  is  the  decay  of  the  explosion  itself.  When  the  afterburn  has  consumed 
all  the  residual  fuel,  the  event  decays  back  to  ambient  conditions.  This  decay  feature 
lasts  between  1-3  seconds.  These  three  features  make  up  the  primary  time  history  of  the 
‘sensed’  explosion.  Figure  1.5  shows  total  radiance  as  a  function  of  time. 
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Figure  1.5  Time  History  of  Typical  Event 

The  three  features  for  the  event  shown  in  Figure  1 .5  are  easily  distinguished. 
Explosion  initiation  occurs  at  1.7  seconds,  followed  by  secondary  afterburn  reaching  a 
peak  value  at  approximately  2.2  seconds  then  the  decay  period  out  to  5  seconds. 

By  assuming  that  a  detonation  can  be  primarily  described  as  a  gray  body,  some  general 
characteristics  of  a  detonation  event  can  be  derived.  The  temperature  of  a  typical 
explosion  initiation  feature  can  be  estimated  as  1800°K  (Cooper,  123).  Similarly,  the 
afterbum  feature  can  be  estimated  as  1 100°K,  and  the  decay  feature  can  be  estimated  as 
400°K.  By  assuming  a  black  body  emitter,  values  for  the  features  of  emitted  spectra  can 
be  predicted.  Using  Equations  1.2,  1.3,  and  1.4,  the  maximum  wavelength,  frequency 
and  the  total  energy  of  the  three  temporal  features  of  a  typical  detonation  event  can  be 
estimated. 
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Table  1.1  Predicted  Features  of  Event  Emitted  Radiation 


Feature 

Temperature 

(Kelvin) 

Wavelength  of 
Maximum 
Intensity  A,m 
(|um) 

Total 

Emittance 

(W/cm2) 

Initiation 

1800 

1.61 

59.50 

Afterburn 

1100 

2.63 

8.30 

Decay 

400 

9.66 

0.05 

Table  1.1  shows  expected  values  for  the  wavelength  of  maximum  emittance  and 
total  emittance  of  each  feature. 

Problem  Statement 

An  objective  of  the  Radiant  Brass  test  program  is  to  explore  the  potential  use  of 
space-based  MASINT  sensors  to  support  combat  units  in  battle-space  characterization. 
Little  infrared  spectral  information  is  available  on  the  signatures  of  conventional 
weapons.  AFIT’s  collection  from  the  Radiant  Brass  tests  exhibited  promise  of 
identification  by  describing  event  signatures.  The  objective  of  this  thesis  will  be  to 
interpret  infrared  spectral  signatures  of  detonation  events  already  collected  and  try  to 
provide  some  insight  into  the  identification  of  ordnance  or  event  conditions.  In  order  to 
accomplish  this  task  a  good  understanding  of  the  collection  process  is  necessary. 
Additionally,  an  understanding  of  the  data  processing  procedures  is  essential.  Finally, 
external  variables  affecting  the  event  signal  collection  must  be  understood.  In  the  case  of 
IR  remote  sensing,  atmospheric  absorption  is  the  largest  factor  affecting  signal  collection 
in  certain  frequencies.  Once  an  atmospheric  adjustment  is  made,  then  the  data  is 
evaluated  at  each  time  interval  in  an  effort  to  find  one  or  more  features  that  can  be  used 
for  identification  of  detonation  type. 
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Scope 

The  focus  of  the  efforts  is  to  determine  the  feasibility  of  using  the  data  collected 
with  an  FTIR  to  discriminate  high  temperature  events.  This  approach  will  lead  to  finding 
one  or  two  features  of  high  temperature  events  that  can  be  used  to  describe  the  event. 

This  thesis  begins  with  the  spectral  and  radiometric  data  from  detonation  events  collected 
during  previous  research.  As  a  feasibility  study,  some  assumptions  are  made  to  reach  the 
end  results.  However,  any  assumptions  will  be  explained  in  detail  before  being-applied. 
Much  of  the  efforts  in  this  research  will  be  on  the  manipulation  of  the  detonation  data  and 
the  PLEXUS  transmittance  data.  The  two  data  types  are  mathematically  merged  and 
evaluated  together.  The  data  was  saved  in  the  American  Standard  Code  for  Information 
Interchange  (ASCII)  format  as  *.txt  files. 

Summary  of  current  knowledge 

Orson  (17)  used  AFIT’s  spectral  radiometer  to  compute  the  total  energy  from 
detonation  events.  Since  there  were  other  agencies  collecting  data  at  the  detonation 
events  with  radiometers,  the  unique  opportunity  to  compare  data  was  available.  Data  was 
integrated  in  four  mostly  atmospherically  transparent  spectral  bands  to  compute  an 
energy  total  for  each  band.  These  four  energy  totals  were  then  compared  to  the  energy 
recorded  in  the  same  spectral  bands  using  radiometers  of  other  agencies  collecting  data. 
This  approach  was  invaluable  in  verifying  the  quality  of  our  equipment,  procedures  and 
data  collected.  This  work  showed  that  the  data  collected  with  our  spectral  radiometer 
was  within  15%  of  the  radiometer  data  of  other  data  collection  teams.  Past  research  also 
discussed  pattern  recognition  options  for  event  discrimination  and  possible  future 
research  options.  This  research  will  build  on  previous  work  and  focus  on  temporal 
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changes  in  event  temperature  and  event  size  for  event  characterization.  Due  to  the 
inherent  difficulties  and  costs  associated  with  bomb  detonation  collection,  little  spectral 
data  is  available  on  conventional  munitions.  The  analysis  in  this  thesis  is  limited  to  the 
spectral  data  collected  from  the  Radiant  Brass  2  A  and  2B  tests  and  signatures  collected 
fr  om  the  Radiant  Brass  3  A  and  3B  tests. 

The  experimental  equipment  used  for  collection  of  the  spectral  signature  data  was 
a  Bomem  FTIR  MR- 154  spectrometer.  This  collection  system  was  comprised  of  two 
primary  components,  the  spectrometer  itself  with  its  associated  optics  and  the  data 
control  computer.  The  spectrometer  was  placed  on  a  tripod  in  clear  view  of  the 
exploding  ordnance.  Optics  on  the  front  of  the  spectroradiometer,  were  used  to  collect 
source  radiation  and  help  control  the  system  field  of  view.  The  field  of  view  of  the 
system  was  approximately  a  300  m  diameter  circle  at  5  km.  Two  detectors,  an  Indium 
Antimonide  (InSb)  and  a  Mercury,  Cadmium,  Telluride  (HgCdTe)  determined  radiance 
information.  A  liquid  nitrogen  dewer  attached  to  the  instrument  provided  a  stable  cold 
reference  source.  A  video  signal  was  monitored  at  the  spectroradiometer  to  improve 
pointing  accuracy.  The  spectroradiometer  was  connected  to  the  data  control  computer  by 
a  serial  cable.  A  feed  from  the  instrument  optics  was  placed  in  the  vicinity  of  the  control 
PC.  This  video  signal  was  used  to  provide  real  time  images  to  the  operator  and  for 
permanent  record  on  Hi-8  digital  tape.  All  measured  signatures  were  stored  on  the 
control  PC  after  each  event.  The  equipment  was  calibrated  on  site  with  a  variable 
temperature  black  body. 

The  Bomem  Acquire  software  processed  the  spectral  signatures.  After  proper 
calibration  and  data  scrubbing  the  spectrum  was  saved  in  ASCII  format. 
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Approach 

The  primary  source  of  data  to  be  used  in  this  research  effort  will  be  the  calibrated 
ASCII  data  from  Radiant  Brass  2  and  Radiant  Brass  3  created  during  previous  research. 
Since  these  are  very  large  data  files  (2200  frequency  bins  by  250  time  bins)  up  to  12 
Megabytes,  FORTRAN  was  used  for  data  manipulation  and  calculations.  Mathematica, 
MicroCal  Origin  and  Table  Curve  2-D  were  used  to  verity  the  results  of  the  FORTRAN 
output  and  as  display  tools. 

Locally  created  atmospheric  transmittance  files  will  be  used  as  a  partial  solution 
to  atmospheric  absorption  of  signal.  PLEXUS  will  be  used  to  ‘undo’  atmospheric 
absorption  as  much  as  possible  to  estimate  the  signal  from  the  source  of  the  emissions. 

In  addition  to  the  continuation  processing  data  collected  during  past  research,  data 
from  various  rocket  launches  were  collected  as  part  of  this  research. 

Summary 

This  completes  the  introduction  material  for  the  collection  and  analysis  of 
conventional  ordnance  spectral  signatures.  Chapter  II  will  provide  a  discussion  of  FTIR 
operations,  data  manipulation,  and  fitting  procedures  used  to  isolate  features  of  the  data. 
Chapter  III  will  detail  the  research  process  step  by  step.  Analysis  performed  are 
presented  in  Chapter  IV.  Data  evaluation  results  are  examined  in  Chapter  V.  Finally, 
Chapter  VI  will  formulate  conclusions  and  offer  recommendations  for  future  study. 
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II.  Theory 


Introduction 

Previous  research  focused  on  the  remote  sensing  collection  and  data  verification 
of  various  data.  This  data  included  detonations  of  various  explosives.  Data  was  collected 
using  AFIT’s  Bomem  FT-154  Fourier  Transform  Infrared  (FTIR)  spectroradiometer. 

The  focus  of  this  thesis  is  to  process  and  interpret  the  IR  emissions  collected  during 
previous  research. 

This  chapter  will  provide  the  basic  theory  needed  for  a  reader  not  familiar  with 
the  data  collection  using  the  Bomem  FT-154  FTIR.  The  first  section  will  provide  an 
understanding  of  basic  remote  sensing  principles.  The  second  section  will  describe  the 
method  of  information  collection,  including  the  FTIR  concept  of  operations,  radiance 
calibration  theory  and  sensor  fusion.  A  review  of  the  sensed  object,  detonations,  will  be 
covered  in  the  bomb  phenomenology  section.  Finally,  an  in  depth  discussion  of  the 
approach  and  specific  methods  used  to  reach  in  my  research  will  be  discussed. 

Since  the  spectroradiometer  software  outputs  all  information  in  cm'1,  these  units 
will  be  used  primarily  in  this  thesis. 

Infrared  Remote  Sensing 

The  region  of  the  electromagnetic  spectrum  used  for  this  work  is  the  infrared 
region.  This  region  of  the  spectrum  is  especially  well  suited  for  detecting  ordnance 
detonations  for  three  reasons.  First,  the  center  frequency  of  spectral  emission  is  centered 
in  IR  region  for  the  afterburn  and  decay  as  mentioned  in  Chapter  I  and  shown  in  Table 
1.1.  Second,  the  IR  region  provides  atmospheric  “windows”  that  allow  the  signal  at 
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certain  frequencies  to  pass  through  the  atmosphere  with  little  or  no  loss.  Finally,  the 
contrast  in  emittance  between  a  detonation  event  and  the  corresponding  IR  background  is 
very  large.  A  high  contrast  allows  for  easier  detection. 

Figure  2.1  shows  the  location  of  the  IR  region  of  the  electromagnetic  spectrum. 
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Figure  2.1  Electromagnetic  Spectrum 

Figure  2.2  shows  the  location  and  type  of  background  emission  sources  as  well  as 
many  commonly  used  IR  terms. 
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Figure  2.2  IR  Spectral  Regions  [pm  (cm1)] 

SWIR  stands  for  Short  Wave  Infrared  while  MWIR  and  LWIR  are  Mid  Wave 
Infrared  and  Long  Wave  Infrared  respectively.  SWIR  is  comprised  of  reflected  solar 
radiation  and  disappears  at  night.  During  the  day,  the  MWIR  is  a  combination  of 
reflected  and  emitted  radiation  while  at  night  it  is  only  emitted.  The  LWIR  is  always 
emitted  radiation. 

By  its  nature  a  detonation  event  is  an  emitter  of  radiation.  The  form  and  intensity 
of  spectral  radiance  from  a  graybody  emitter  can  be  produced  using  Equation  1.1.  Since 
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the  amount  of  energy  released  in  a  detonation  is  large,  there’s  a  sharp  contrast  to  the 
relatively  low  emittance  from  the  ambient  background. 


In  the  discussion  of  electromagnetic  energy  emission,  propagation  and  detection, 
a  number  of  quantities  are  commonly  used.  Table  2.1  defines  a  number  of  radiometric 
quantities  used  to  describe  the  energy  content  of  incoherent  radiation  fields. 


Table  2.1  Radiometric  Quantities  and  Units  (21, 14-26) 


Quantity 

Symbol 

Units 

Equation 

Defining  Equation 

Radiant  energy 

Q 

Joules 

2.1 

O 

ii 

e 

& 

Radiant  energy  density 

w 

Joule/cm3 

2.2 

wJS. 

dV 

Radiant  flux 

0) 

watt 

2.3 

dt 

Radiant  flux  density 
(Irradiance) 

E 

watt/m2 

2.4 

r 

E  = - 

dt 

Radiant  intensity 

I 

watt/str 

2.5 

j  d® 

dt 

Radiance 

L 

watt/cm2  str 

2.6 

L-?l 

dA 

A  more  detailed  description  of  the  spectral  and  the  energy  measurement  methods 
performed  by  the  FTIR  spectroradiometer  will  be  described  in  the  next  section. 

Fourier  Transform  InfraRed  (FTIR)  Concept  of  Operation 

The  detonation  data  used  during  this  thesis  was  taken  using  a  Bomem,  model  MR- 
154,  Fourier  Transform  Infrared  spectroradiometer.  Since  all  data  was  collected  and 
processed  by  this  instrument,  it  is  important  to  describe  the  operation  of  the  equipment  in 
detail. 
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Many  methods  for  obtaining  spectral  information  rely  on  the  dispersion  of  light  to 
achieve  spectral  separation.  This  is  commonly  done  with  a  prism  or  grating.  The  MR- 
154  uses  a  Michelson  interferometer  to  create  an  interferogram  of  the  input  signal.  A 
fourier  transform  of  the  interferogram  creates  the  spectra.  One  advantage  of  the  FTIR 
technique  is  all  frequencies  of  the  input  source  are  falling  on  the  detector  at  all  times. 
Figure  2.3  shows  the  path  of  the  collected  IR  radiation  inside  the  interferometer. 


Figure  2.3  MR-154  Scan  Arm.  (5, 2-5) 

The  input  beam  in  Figure  2.3  is  gathered  by  the  collection  optics  and  collimated 
by  the  optics  which  focus  the  input  beam.  Upon  striking  the  beam  splitter,  the  incident 
beam  is  split  into  two  different  paths.  As  the  scan  arm  moves,  the  optical  path  length 
changes  between  the  two  beams.  The  two  beams  are  recombined  to  create  constructive 
and  destructive  interference  patterns  based  on  the  differences  in  the  optical  path  lengths 
Generally,  intensities  for  all  arm  positions  in  both  scan  directions  are  added  together  to 
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get  one  interferogram.  These  intensities  are  subsequently  decomposed  into  a  spectrum  by 
Fourier  transforming  the  interferogram.  In  normal  operations,  one  scan  of  the  Michelson 
interferometer  is  defined  as  the  forward  and  reverse  sweep  of  the  scan  arm  In  zero  scan 
mode,  a  scan  is  only  one  sweep  direction. 

The  nature  of  the  FTIR  collection  method  provides  inherent  broadband  capability. 
No  optical  filtering  is  used  and  the  spectral  bandwidth  is  dependant  on  the  detector  and 
beam  splitter  not  the  system  itself.  This  provides  great  flexibility  as  different  detectors 
can  be  used  with  different  frequency  responses.  The  maximum  optical  path  length 
difference  of  the  system  limits  the  spectral  resolution.  Within  the  system  limits,  the 
control  computer  can  make  quick  resolution  changes. 

Digitizing  of  the  interferogram  requires  precise  monitoring  of  the  optical  path  in 
the  interferometer.  The  MR-154  uses  an  internal  He-Ne  laser,  15798  cm'1,  to  determine 
the  optical  path  difference  feedback  for  the  system  (5,  2-6).  The  laser  provides  a  known 
reference  signal.  This  reference  ensures  very  stable  optical  path  length  measurements.  A 
liquid  nitrogen  cold  reference  source  is  attached  to  the  MR- 154  for  an  accurate  cold  body 
reference.  Determining  absolute  intensity  information  is  a  function  of  the  cold  reference 
source,  choice  of  detector,  and  most  importantly  calibration. 

The  MR- 154  collects  radiant  energy  as  a  function  of  frequency.  Spectral  radiance 
can  be  derived  from  the  signal  when  the  detector  hardware  and  software  are  given  good 
calibrations  and  input  parameters.  A  properly  set  up  and  well-calibrated  system  provides 
useful,  accurate  measurements  that  represent  the  emissions  from  an  unknown  source. 

The  Bomem  Acquire  software  controls  the  Bomem  MR- 154.  During  operation, 
one  interferogram  is  collected  using  the  forward  and  reverse  scans  of  the  Michelson 
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interferometer.  Each  measurement  is  a  number  of  interferograms  averaged  together.  The 
number  of  interferograms  to  average  is  selected  by  the  user.  Acquire  also  controls  the 
spectral  range  of  each  detector,  the  number  of  measurements  to  be  collected,  and  the  time 
interval  between  measurements. 

The  process  of  creating  an  accurate  measurement  of  an  unknown  has  two  steps, 
calibration  and  unknown  acquisition.  The  first  step  is  to  create  calibration  files  using 
known  reference  sources.  A  good  calibration  consists  of  at  least  two  radiometric 
reference  points.  Ideally,  the  reference  points  should  straddle  the  anticipated  unknown 
temperature.  Each  reference  point  is  built  with  the  same  user-selected  variables  used  for 
measurement.  A  radiometric  reference  file  is  created  for  each  reference  point.  This 
process  is  done  by  the  Bomem  Acquire  software  package  used  with  the  instrument. 

The  calibration  and  data  collection  process  is  detailed  in  Figure  2.4 
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Figure  2.4  Software  Measurement  Process  (5,  7) 

The  second  step  to  creating  accurate  measurements  from  an  unknown  source  is 
data  acquisition.  This  is  accomplished  directly  in  the  Acquire  software  program.  The 
user  selects  spectral  range,  number  of  scans  to  average  per  measurement,  and  number  of 
measurements.  The  system  collects  the  event  data  and  Acquire  creates  an  interferogram 
file.  To  process  the  interferogram,  a  cosine  apodization  function  is  applied.  This 
corrected  interferogram  is  transformed  via  a  Fast  Fourier  Transform  (FFT),  to  create  a 
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complex  raw  spectrum  file.  Applying  the  radiometric  calibration  file  to  the  complex  raw 
spectrum  file  creates  the  final  spectral  radiance  measurement. 

Each  block  in  Figure  2.4  containing  a  bracketed  quantity  represents  a  specific 
output  file  from  the  Acquire  software.  The  bracketed  quantities  refer  to  the  file  extension 
attached  to  a  user-defined  filename.  An  ‘x’  refers  to  the  detector  position,  A  or  B,  which 
is  being  processed.  The  Acquire  software  also  provides  limited  data  manipulation 

Atmospheric  effects 

Any  time  electromagnetic  energy  propagates  through  the  atmosphere,  it  is 
affected  by  atmospheric  absorption  and  scattering.  The  atmosphere  is  comprised  of  many 
different  molecules.  Each  molecule  will  absorb  different  frequencies  of  electromagnetic 
radiation  due  to  quantum  mechanical  interactions.  These  quantum  mechanical 
interactions  include  electronic,  vibrational  and  rotational  modes  of  the  molecule.  Due  to 
relatively  low  energy  levels  associated  with  IR  radiation,  electronic  transitions  are  not 
prevalent.  The  vibrational  and  rotational  modes  are  the  strong  absorbing  mechanisms. 
The  molecules  that  absorb  electromagnetic  radiation  in  the  IR  spectrum  include  H2O, 

CO2,  and  O3.  Figure  2.5  summarizes  the  typical  vertical  atmospheric  absorption  profile 
for  the  Earth  atmosphere. 
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H20  and  C02  have  many  absorption  bands.  Each  large  band  is  related  to  a  vibrational 
mode.  The  width  of  the  absorption  is  due  to  the  rotational  modes.  The  03  absorption 
spectrum  around  10  pm  is  not  an  issue  for  this  thesis  since  03  is  a  minor  gas  below  100 
km  in  altitude. 

Concentrations  of  H20,  C02,  and  03  are  dependent  on  atmospheric  conditions. 
These  concentrations  vary  from  day  to  day.  The  amount  of  absorption  is  a  function  of 
atmospheric  composition  and  the  distance  to  the  source.  Atmospheric  absorption  is 
commonly  described  with  a  transmission  coefficient  between  zero  and  one.  A  coefficient 
of  one  means  no  attenuation  and  coefficients  near  zero  are  totally  attenuated. 

The  total  attenuation  of  the  atmosphere  is  a  combination  of  absorption  and 
scattering.  Particulates  in  the  atmosphere  with  radii  between  0. 1  and  10  pm  attenuate 
signal  between  visual  and  thermal  IR  frequencies  (9,  279).  The  amount  of  scattering  is 
based  on  the  particulates  cross-sectional  area  and  the  frequency  of  the  source  radiation. 
Examples  of  these  particle  sources  are  clouds,  smoke,  and  dust. 
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Data  Manipulation 

Data  files  created  from  detonation  events  are  very  large.  Depending  on  spectral 
resolution  and  data  collection  time,  a  single  file  may  have  as  many  as  300  temporal  steps 
with  up  to  2500  frequency  steps.  FORTRAN  was  selected  for  its  ability  to  handle  the 
large  data  files  and  manipulate  the  data  efficiently.  FORTRAN’S  strength  is  in  building 
the  data  sets  into  two-dimensional  arrays.  A  2-D  array  allows  interrogation  of  the  data  by 
timestep  or  frequency.  In  building  the  code,  great  care  was  taken  to  allow  for  different 
sizes  of  data  files.  It  was  important  to  have  a  program  that  could  determine  the  size  of 
the  data  file  and  read  the  data  regardless  of  the  size  or  shape  of  the  ASCII  file  being 
accessed.  The  program  had  to  read  the  ASCII  files  created  from  the  detonation  data  and 
the  PLEXUS  transmittance  profiles.  The  FORTRAN  code  was  modified  to  read 
PLEXUS  files  and  interpolate  the  value  of  the  transmittance  values  for  the  same 
frequencies  as  the  ASCII  data. 

Once  the  interpolation  was  complete,  the  ASCII  data  could  be  divided  by  the 
transmittance  value  at  each  frequency  to  attempt  recreation  of  the  signal  amplitude  at  the 
source.  If  the  PLEXUS  value  at  a  given  frequency  were  too  small,  the  adjusted 
amplitude  value  at  that  frequency  would  be  set  to  zero.  Setting  the  value  for  adjusted 
intensity  to  zero  when  the  transmittance  was  below  a  cutoff  value  had  two  advantages. 
First,  setting  the  adjusted  value  at  that  frequency  to  zero  avoided  overestimated  values 
when  dividing  by  a  very  small  number.  Second,  if  the  transmittance  is  small,  the 
confidence  in  any  signal  collected  at  that  frequency  would  also  be  small.  This  method 
allowed  any  points  with  adjusted  values  of  zero  to  be  ignored  when  fitting  to  a  function. 
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Mathematica  was  linked  to  the  FORTRAN  program  to  allow  for  graphical 
interface.  Mathematica  allowed  the  user  to  edit  the  location  and  filenames  of  the  input 
ASCII  intensity  data  and  PLEXUS  transmittance  files  and  the  output  data  files  created  by 
FORTRAN.  After  the  program  was  run,  Mathematica  could  quickly  create  graphical 
representations  of  the  data  created.  The  union  of  Mathematica  and  FORTRAN  proved  to 
be  an  effective  combination  of  their  respective  graphical  and  computational  strengths. 
With  the  program  reading  the  input  files,  calculating  new  values  and  creating  adjusted 
data  files  effectively,  the  identification  of  temporal  features  of  an  event  could  finally 
begin. 

Data  Evaluation 

Adjusted  intensity  data  would  need  to  be  evaluated  for  each  time  step  within  an 
event  data  file.  The  goal  of  this  evaluation  was  to  find  one  or  two  meaningful  parameters 
for  each  time  step  that  could  be  followed  over  the  course  of  the  data  collection.  The 
evolution  of  the  selected  parameters  as  a  function  of  time  were  then  investigated  for 
identifiable  features  unique  to  each  type  of  event.  The  approach  selected  was  to  fit  event 
temperature  and  size  to  the  Planckian  curve  representation  of  two  blackbodies.  Event 
size  and  temperature  were  adjusted  during  fitting. 

Chapter  III  will  detail  the  derivation  of  the  equations  used  to  find  event 
temperature  and  event  size.  However,  an  overview  of  their  derivation  is  necessary  prior 
to  a  discussion  of  curve  fitting  methods  used. 

The  general  approach  used  here  will  be  to  assume  the  signal  recorded  by  the 
instrument  is  a  result  of  a  single  gray  body  emitter  or  the  combination  of  two  blackbody 
emitters.  The  Planck  function  relates  temperature  and  frequency  to  the  radiation  of  a 
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gray  body.  This  function  was  introduced  in  chapter  I  as  Equation  1.5  Chapter  III  will 
detail  the  assumption  of  a  black  body.  The  substitution  of  e(v)  =  1  creates  the  new  form 
Equation  2.1 


cv  t\  2hc  v  /"o  i 

S(vJ)=-7iZ7 -  v2'1 

e  /kT-\ 

Prior  to  detonation,  the  detector  is  collecting  data  from  the  scene  within  the  field 


of  view.  This  will  be  called  the  background  signal  and  named  Sback  with  a  temperature 


Tback.  This  signal  fills  the  entire  field  of  view.  Since  the  atmospheric  temperature  at 
detonation  time  is  a  known  quantity,  the  form  of  the  signal  from  the  background  is  known 
for  all  frequencies.  Calibration  measurements  were  created  using  a  blackbody  that  didn’t 
completely  fill  the  field  of  view.  The  field  of  view  filled  by  the  calibration  blackbody 
relative  to  the  entire  field  of  view  occupied  by  Shack  changes  the  amplitude  of  the  signal. 

A  new  variable,  K,  is  introduced  to  represent  this  number.  The  other  blackbody  source  of 
radiation  is  the  detonation  event.  The  signal  from  this  source  will  be  named  Sbomb  with  a 
temperature  Tbomb .  The  bomb  occupies  a  dynamic  fraction  of  the  field  of  view  during 
initiation,  afterburn  and  decay.  The  fraction  of  the  field  of  view  occupied  by  the 
detonation  will  be  represented  by  the  new  variable,  F.  The  total  signal  in  the  field  of 
view  of  the  detector  is  then  called  SSCene .  All  terms  are  combined  into  Equation  2.2  to 


describe  the  signal  from  the  scene  at  a  given  time  step. 


SScene=K[F-Sbomb  +  (l-F)-Sbb\  (2.2) 


This  is  the  basic  form  of  the  Equation  that  will  be  used  to  describe  the  detonation 
event  recorded  by  the  instrument.  All  quantities  are  known  for  each  time  step  except 
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fractional  field  of  view  filled  by  the  detonation  and  bomb  temperature.  These  become  the 
two  variables  found  during  the  curve  fitting  process. 


Fitting  Procedures 


Since  Chapter  III  will  provide  a  detailed  review  of  the  fitting  procedures  used, 
only  an  overview  of  the  process  will  be  discussed  here.  Fitting  two  variables 
simultaneously  in  an  exponential  Equation  proved  to  be  challenging.  A  modified  method 
of  bisection  was  created  to  fit  the  two  variables.  Event  size  was  selected  as  the  first 
variable  evaluated.  Root  Mean  Square  Error  (RMSE)  was  used  as  the  convergence 
criteria.  RMSE  is  described  as  the  square  root  of  the  sum  of  residuals  squared  divided  by 
the  number  of  residuals.  The  calculation  of  RMSE  for  a  single  timestep  is  shown  in 
Equation  2.18 


RMSE  = 


f  nw  a 

Xtv,.(z,.-z,.)2 


(=1 


\ 

/nw-  2 

J 


(2.18) 


where  i  is  the  frequency  location,  w,-  is  the  weighting  of  each  data  point,  z(.  is  the  data 
estimate  using  variables,  z,  is  the  actual  data  being  fit,  and  nw  is  the  total  number  of  data 
points  used.  The  term  nw-2  is  the  Degrees  Of  Freedom  (DOF).  For  this  research, 
weighting  will  be  accomplished  during  PLEXUS  data  manipulation  so  w,  is  1  for  all 
wavenumbers.  Equation  2.18  is  minimized  during  curve  fitting  using  the  bisection 
method  described  below. 

The  method  of  bisection  fits  an  Equation  by  changing  variables  in  decreasing 
intervals.  The  method  uses  a  repeated  halving  of  subintervals.  The  half  with  the  smallest 
error  is  then  identified.  Equation  2. 19  shows  the  interval  magnitude  for  variable  m  at 
iteration  n. 
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Interval ”  =  (0.5)"  •  Varguessm 


(for  n  >  1) 


(2.19) 


For  each  interval,  the  RMSE  for  Varguess,  as  well  as  the  values  of  Varguess  + 
interval  are  found.  The  value  with  the  lowest  RMSE  becomes  the  new  Varguess.  This 
method  is  repeated  until  an  established  criteria  is  met.  Many  methods  in  print  use  a 
maximum  number  of  iterations  as  the  exit  criteria.  To  ensure  the  most  accurate  variable 
values,  bisection  was  continued  until  three  consecutive  iterations  with  both  variables 
yielded  no  change  in  RMSE.  Other  statistical  values  were  computed  as  well.  These 
values  were  computed  to  compare  the  fit  of  the  FORTRAN  code  to  the  fit  of  the  data 
using  TableCurve.  These  values  and  their  derivation  are  described  in  detail  next  chapter. 
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III.  Approach 


Introduction 

Previous  research  provides  a  source  of  well-conditioned,  calibrated  data  from 
explosive  detonation  tests.  Through  radiometric  comparisons  to  data  of  the  same  event 
collected  by  other  agencies,  the  validity  of  the  data  was  verified.  The  approach  taken  in 
this  research  will  be  to  use  the  existing  data  to  find  temporal  profiles  of  event  size  and 
temperature  for  each  event.  Since  data  taken  from  the  static  detonations  of  Radiant  Brass 
3B  are  the  most  consistent,  this  data  set  is  used. 

Radiant  Brass  3B  consisted  of  23  statically  detonated  events  over  the  course  of 
two  nights.  Three  different  explosive  types  and  three  different  explosive  amounts  were 
used  during  this  test.  The  ordnance  was  situated  on  solid  ground  and  propped  up  on 
wooden  tripods  at  a  45  degree  elevation  angle  tail  high.  Each  bomb  face  was  pointed 
toward  the  instrumentation.  The  larger  bombs  were  placed  in  clay  craters  tail  high  at  the 
largest  elevations  possible 

This  chapter  begins  with  a  discussion  of  the  data  collection  and  calibration 
procedures  used  to  create  the  data  used  in  this  research.  A  more  detailed  discussion  on 
the  collection  and  calibration  of  this  data  can  be  found  in  the  BOMEM  Users  Manual,  (2 
and  3)  or  in  Orson’s  thesis  (17).  This  chapter  will  then  explain  the  integration  of  the 
programs  Mathematica,  PLEXUS,  FORTRAN  and  Table  Curve  used  to  process  and 
interpret  the  data.  The  remainder  of  chapter  will  discuss  the  development  and 
implementation  of  the  numerical  methods  used  to  fit  the  data. 
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Calibration 


During  data  collection,  great  care  was  taken  to  obtain  the  best  data  possible. 
Accurate  calibration  references  were  created  using  a  calibrated  black  body  with  a  known 
aperture  setting.  The  blackbody  was  positioned  a  distance  from  the  detector  that 
simulated  the  size  of  the  event  in  terms  of  fractional  field  of  view.  To  avoid  detector 
saturation,  the  anticipated  fractional  field  of  view  seen  by  the  detector  at  bomb  initiation 
was  used.  From  previous  data,  it  was  determined  the  size  of  a  detonation  event  at  4750  m 
could  be  approximated  by  a  0.2”  aperture  at  76”.  The  blackbody  was  set  at  it’s  maximum 
value  of  1000°C.  In  this  configuration,  the  MR- 154  focused  on  the  black  body  to 
determine  the  optimal  gain  settings.  The  detector  gain  settings  were  selected  to  produce 
maximum  signal  available  without  detector  saturation.  The  final  configuration  settings 
weie  selected  to  maximize  spectral  range  of  the  collected  event  signature  with  both 
detectors.  The  spectral  range  of  the  Mercury  Cadmium  Telluride  (HgCdTe)  detector  was 
set  at  500  to  6000  cm1.  The  spectral  range  of  the  Indium  Antimonide  (InSb)  was  set  at 
1800  to  6000  cm'1.  These  spectral  ranges  cover  most  of  the  SWIR,  the  MWIR  and  LWIR 
regions.  In  this  collection  configuration  both  detectors  have  some  overlap  in  spectral 
range  for  comparison  between  detectors.  All  event  measurements  were  taken  in  zero 
scan  mode  utilizing  only  one  sweep  direction  of  the  Michelson  Interferometer.  Finally, 
events  were  to  be  collected  at  16  cm 1  spectral  resolution  with  temporal  resolution  of 
0.047  seconds  and  at  4  cm'1  spectral  resolution  and  with  resolution  of  0.123  seconds. 

Once  gain  settings,  temporal  resolution,  spectral  resolution  and  data  collection 
times  were  set,  calibration  files  were  created.  Calibration  files  were  made  with 
radiometric  reference  files.  A  calibrated  blackbody  source  was  used  to  build  the 
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radiometric  reference  files  at  700°C  and  900°C.  A  Bomem  Optics  74  milliradian 
collimator  optics  telescope  was  used  during  Radiant  Brass  3B  data  collection.  This 
telescope  resulted  in  a  roughly  360  m  diameter  field  of  view  at  the  target  point. 

The  Bomem  MR  Series  FT-Spectroradio meter  Design  Overview  and  Theory 
section  of  the  Bomem  Users  Manual  (4)  details  the  calibration  and  data  collection  process 
of  the  MR- 154.  The  appropriate  measurement  equation  for  a  blackbody  source  is  shown 
as  Equation  3.1.  The  units  of  Smeas(v)  is  W/cm2  str  cm'1  and  will  be  called  scene  spectral 
radiance. 

SmeJv)  =  AoptAhh  ^f-[L(v)Det(v  )Tropt  (v  )Trfilter  (y)Tram  (v )  ]  (3.1) 


where 


Aopt  =  area  of  instrument  optical  aperture 

Abb  =  area  of  blackbody  emitter 

R  =  distance  from  calibration  source  to  instrument 

L(v)  =  blackbody  spectral  radiance  W/cm2  sr  cm'1 

Det(v)  =  detector  frequency  dependent  power  responsivity 

Tropt(v)  =  instrument  optics  wavelength  dependent  transmission 

Trfiiter(v)  =  spectral  filter  transmission 

Tratm(v)  =  atmospheric  transmission 

Geiec  =  electronic  gain 

The  MR- 154  has  no  filter  so  Trfliter(v)  can  be  neglected.  All  system  wavenumber 
dependant  terms,  Det(v)  and  Tropt(v),  as  well  as  Aopt  and  Geiec  can  be  grouped  into  one  on 
term K(v).  With  these  manipulations  Equation  3.1  becomes  Equation  3.2: 

S„.Jv)  =  |r*(v)[. I(v)7rra(v)]  (3.2) 

Further  assuming  the  transmission  losses  due  to  the  atmosphere  are  negligible  due 
to  the  short  distance  eliminates  Tratw(v).  For  a  calibration  sequence  the  Abt/R2  is  constant 
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and  can  be  combined  into  K(v).  Lastly,  when  the  event  or  calibration  does  not 
completely  fill  the  system  field  of  view,  L(v)  can  be  expanded  to: 

Su,Jv)  =  *(v)[Z*(v)  +  (v)  +  Ms"»]  (3.3) 


where 

Lbb(v)  =  black  body  source  radiance 

Lback(v)  =  source  background  radiance 

MStray( v)  =  system  stray  radiance  term 

To  eliminate  the  source  background,  two  measurements  are  taken  at  each 
calibration  reference  temperature,  a  source  plus  background  and  a  background.  These 
two  reference  measurements  are  subtracted  to  eliminate  Lback(v).  Theoretically,  this 
leaves  only  the  contribution  of  the  black  body  (bb)  source  in  the  S,  L,  and  M  quantities. 
Once  Lback(v)  is  subtracted,  the  equation  becomes: 

Sbb(v)  =  K(v)(Lbb(v)  +  MStray(v))  (3.4) 

A  proper  calibration  procedure  requires  collection  of  a  minimum  of  two  reference 
temperatures,  commonly  called  hot  (H)  and  cold  (C).  These  reference  temperature  files 
each  need  a  separate  background  reference.  The  background  reference  theoretically 
removes  all  signal  not  associated  with  the  calibration  source.  To  create  the  background 
references,  the  blackbody  radiation  source  is  removed  from  the  field  of  view.  Any 
radiation  in  the  field  of  view  is  captured  and  then  subtracted  from  the  hot  and  cold 
reference  files.  With  these  two  background  subtracted  reference  measurements,  Aftray( v  ) 
and  K(v )  are  found  by  Equations  3.5  and  3.6. 


^(v)=Mv)-*c(v) 

Lh(v)-Lc(v) 


(3.5) 
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1/ Stray  /  \  _  L H  (y)Sc(y)-Lc(y)SH(y) 
SsH(v)~S?eas(v ) 


(3.6) 


Once  Aftray(v )  and  K(v  )  are  determined,  the  equation  to  calibrate  successive 
measurements  is: 

Scdibrated(v)  =  SMeas  (v)r‘(v)-M%(v)  (3.7) 

Equation  3.7  is  essentially  a  mathematic  correction  applied  to  a  ‘sensed’  unknown  signal. 
This  assumption  used  during  calibration  is  only  valid  if  all  operating  parameters  of  the 
system  remain  constant  between  reference  point  acquisition  and  the  unknown 
measurement.  Since  the  MR-154  collects  data  on  both  sweep  directions  of  the 
interferometer,  two  sets  of  calibration  constants  are  created. 

Linearity  is  a  big  assumption  in  Equation  3.7.  This  assumption  requires  linearity 
in  the  system  and  detector.  For  the  large  radiance  and  spectral  range  required  to  measure 
a  detonation  event,  linearity  in  the  InSb  detector  is  a  valid  assumption.  For  the  HgCdTe 
detector  a  quadric  calibration  is  more  appropriate.  The  details  of  a  quadratic  calibration 
can  be  found  in  the  Bomem  users  manual  (5,  3: 109-1 1 1). 

Data  Collection 

The  collection  of  data  of  an  unknown  event  is  considered  in  much  the  same  way 
as  the  radiometric  reference  was  taken  during  calibration.  The  unknown  event  is 
recorded  using  the  same  gain,  resolution  and  detector  settings  used  during  calibration. 

All  of  the  assumptions  that  were  made  between  Equation  3.1  and  Equation  3.2  are  still 
valid.  However,  the  measurement  of  a  dynamic  unknown  blackbody  source  at  a  great 
distance  changes  the  Equation  beyond  the  form  shown  in  Equation  3.2.  Equation  3.2  is 
recalled  as  Equation  3.8. 
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(3.8) 


5,.„,(v)  =  4t/f(v)[L(v)7-r„„(v)] 

For  measurement  of  an  unknown  signal,  all  of  the  radiation  the  detector  is 
receiving  come  from  within  the  field  of  view.  The  detector  signal  received  will  become 
scene  spectral  radiance  represented  by  the  term  SSCene(v  ).  The  term  Tram( v  )  becomes 
very  important  as  distances  increase  and  can  not  be  discarded  as  it  was  during  calibration. 
Chapter  II  discussed  the  effects  of  atmospheric  absorption  across  an  IR  region  of  the 
spectrum.  This  term  remains  in  the  Equation.  The  atmospheric  absorption  will  be 
removed  as  much  as  possible  using  PLEXUS.  The  term  Abt/R2  remains  constant  as  the 
portion  of  the  field  of  view  considered  in  the  radiometric  reference  files.  This  quantity 
can  be  defined  as  the  solid  angle  between  the  detector  and  the  blackbody  aperture  size. 
This  number  can  be  calculated  and  will  be  used  later.  For  now,  it’s  considered  a  constant 
and  is  included  in  K(v).  L(v  )  is  expanded  into  the  radiation  sources  expected  within  the 
FIELD  OF  VIEW  of  the  detector.  For  a  detonation  event,  these  terms  are  radiation  from 
the  explosion,  Sbomb(v ),  and  radiation  from  everything  else  in  the  FOV,  Sback(v  )■  Prior  to 
bomb  initiation,  the  portion  of  the  field  of  view  filled  by  detonation  signal  is  zero.  After 
detonation  the  fraction  becomes  a  number  that  changes  with  time  as  the  event  goes 
through  the  phases  of  initiation,  afterbum  and  decay.  To  account  for  the  dynamic  portion 
of  the  field  of  view  filled  by  the  detonation  event,  a  new  variable  F  is  introduced.  F 
represents  the  relative  fraction  of  the  field  of  view  occupied  by  the  explosion.  The  rest  of 
the  field  of  view  is  then  represented  by  (1-F).  The  new  form  of  the  function  now 
becomes  Equation  3.9. 

^,ie(v)  =  Tratm(v)-X(v)[F 

^ 'bomb  (v)  +  (l -F)  ^ 1  back  (v)]  (3.9) 
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This  Equation  represents  the  signal  received  by  the  detector.  This  total  is  the  sum 
of  the  bomb  radiation  and  background  black  body  radiation.  The  term  SSCene(o)  represents 
the  calibrated  spectral  scene  radiance  data  created  by  the  instrument  Aquire  software.  To 
extract  detonation  temperature  and  size  at  each  time  step,  the  wavenumber  dependant 
Planck  distribution  introduced  as  Equation  1.5  is  recalled  as  Equation  3.10.  Equation 
3.10  describes  the  intensity  and  form  of  spectral  radiance  from  a  gray  body  as  a  function 
of  wavenumber  and  temperature. 

2/zcV3 

S(v,r)  =  e(v)-jc7 -  (3.10) 

e  /kT-l 

The  terms  Sback(v )  and  Sbomb(v )  in  Equation  3.9  are  replaced  with  their 
respective  blackbody  forms  of  Equation  3.10  Gray  body  effects  are  placed  in  the  term 
K(v )  removing  the  e(v )  term  from  Equation  3.10  The  temperature  used  in  Sback( v)  is 
represented  by  T back  and  will  be  the  ambient  atmospheric  temperature  at  data  collection 
time.  The  temperature  in  Sbomb(v )  is  given  the  name  Tbomb  and  is  one  of  the  two  numbers 
to  be  found. 

Finally,  an  atmospheric  transmittance  profile  is  created  using  PLEXUS.  The 
transmittance  profile  uses  the  atmospheric  temperature,  atmospheric  pressure,  humidity 
and  distance  to  detonation  event  to  create  a  file  representing  atmospheric  transmittance  at 
each  wavenumber  location.  If  the  transmittance  value  was  below  0.4,  this  division 
produced  excessive  scene  spectral  radiance  values.  All  wavenumber  locations  with  a 
transmittance  value  of  0.4  or  less  were  considered  less  accurate  than  other  frequencies 
and  scene  spectral  radiance  was  set  to  zero  for  that  frequency.  New  scene  spectral 
radiance  values  were  created  at  wavenumber  locations  having  transmittance  values 
greater  than  0.4  by  dividing  by  the  transmittance  value.  The  PLEXUS  corrected  data  is 
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identified  by  the  *  superscript.  The  form  representing  transmittance  corrected  spectral 
scene  radiance  at  one  time  step  becomes  Equation  3.1 1: 


.f 


2hcv 


2,  .3 


2  hcV 


S(vtene  =  K(v)  ■  F  —  -  -  "  +  K(y)  •  (1  -  F)  . t . t  (3 ■  1 1) 

[Exp{chv  !kTbomb)-  1J  [£jp(cM'/*rtac4)-lJ 


The  asterisk  on  S(v)*scene  indicates  the  original  data  has  been  PLEXUS  corrected. 
This  equation  represents  the  data  for  each  time  step  evaluated  over  all  frequencies.  The 
calibration  term  K(v )  is  distributed  so  that  each  term  may  be  discussed  independently. 

Consider  the  first  term  of  Equation  3.11.  This  term  represents  the  radiation 
received  by  the  detector  due  to  the  detonation  event.  Prior  to  detonation,  F  equals  zero 
and  the  only  radiation  source  in  the  field  of  view  is  the  background.  At  any  time  after 
bomb  initiation,  F  represents  the  fractional  field  of  view  of  the  event  at  that  time  step. 
This  value  is  relative  to  the  portion  of  the  field  of  view  filled  by  the  black  body 
references  during  calibration.  By  relating  the  value  of  Fto  the  known  fraction  of  the 
field  of  view  filled  by  the  reference  blackbody,  event  size  can  be  determined 

When  the  temperature  of  the  event  falls  between  or  near  the  hot  and  cold 
radiometric  reference  files,  the  value  of  K  is  well  described.  Recall  for  all  calibrations 
during  RB3B,  hot  and  cold  references  were  created  at  1173  and  973  Kelvin  respectively. 

The  second  term  in  Equation  3.1 1  describes  the  scene  radiance  values  created  by 
the  Acquire  software  due  to  background  radiation.  Before  detonation,  F  equals  zero  and 
all  radiation  received  is  due  to  the  background  scene  for  that  time  step.  Recall  that  K  is  a 
linear  calibration  parameter  created  by  recording  radiation  from  known  sources  and 
subtracting  background  radiation.  Therefore,  a  perfect  calibration  over  all  temperatures 
would  yield  a  zero  value  for  all  frequencies  when  there  is  no  signal  other  than 
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background.  However,  calibrated  scene  radiance  values  observed  are  not  zero.  For 
radiation  due  to  background,  K  may  not  be  an  appropriate  value  since  the  ambient 
temperature  of  around  300  Kelvin  is  not  between  or  even  near  the  hot  and  cold 
radiometric  reference  temperatures  of  1 173  and  973  Kelvin.  In  order  to  correct  the 
spectral  scene  radiance  values  for  all  data,  average  values  of  the  background  prior  to 
bomb  initiation  are  computed.  Since  the  signal  from  the  event  is  two-three  orders  of 
magnitude  smaller  than  the  signal  of  the  background,  this  can  be  done  without  adversely 
affecting  detonation  data.  The  background  average  values  are  subtracted  from  the 
PLEXUS  adjusted  data  to  create  a  new  data  set  representing  scene  radiance.  With  this 
correction,  the  radiation  at  the  detector  is  described  by  the  new  Equation  3.12. 


S(v,F,Tbomby*ce„e=F< 


2 he  v 


2,  ,3 


Exp(chv  /  kTbomb)  - 1 


(3.12) 


The  term  S(v,F,Tbomb)**cene  has  two  asterisks  representing  the  two  corrections 

made  to  the  data,  PLEXUS  transmittance  correction  and  background  signal  subtraction. 
For  a  given  time  step,  the  equation  has  only  two  variables  when  considered  over  all 
frequencies.  This  is  the  final  form  of  the  function  used  to  find  bomb  temperature  and 
detonation  size  at  each  time  step. 

Data  Interpretation 

Once  corrections  were  made  to  the  data  removing  atmospheric  absorption  and 
background  signal,  the  corrected  data  was  used  to  describe  a  detonation  event.  The  two 
variables  in  Equation  3.12  are  bomb  temperature  and  fractional  field  of  view.  The  data 
was  fit  using  a  simple  least-squares  minimization  procedure.  Curve  fitting  two  variable 
in  a  non-linear  equation  using  least  squares  can  be  accomplished  a  number  of  ways.  The 
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method  chosen  in  this  work  was  bisection  for  two  variables  consecutively.  Non-linear 
fitting  is  an  iterative  procedure  that  begins  with  an  initial  set  of  estimates  of  the  variables. 
Regardless  of  the  method  used  to  converge  to  the  minimum  least-squares  solution,  a 
point-by-point  sum  of  squared  residuals  must  be  calculated  for  each  iteration.  Recall  that 
the  data  correction  method  described  previously  included  subtraction  of  the  average 
background  signal  from  all  time  steps.  Therefore,  any  time  step  with  no  event  in  the  field 
of  view  essentially  contained  no  data.  Curve  fitting  procedures  were  not  applied  to  the 
data  for  time  steps  prior  to  detonation.  At  bomb  initiation,  an  initial  estimate  for  bomb 
temperature  of  1500  Kelvin  was  used.  The  initial  estimate  for  relative  fractional  field  of 
view  was  set  as  the  unit  less  factor  2.  Each  iteration  of  the  curve-fitting  procedure  would 
then  calculate  the  least-square  error  for  the  function  using  the  initial  values  stated  above. 
The  form  of  the  least-square  error  used  for  convergence  is  the  Root  Mean  Square  Error 
(RMSE).  The  form  of  the  RMSE  is  shown  in  Equation  3.13. 


RMSE  = 


(3.13) 


where  j  represents  the  frequency  location  of  the  corrected  data.  The  number  of  frequency 
locations  containing  data  is  nw.  s*J  represents  the  value  of  the  corrected  data  at  each 


frequency.  The  denominator  (nw-2)  represents  the  Degrees  of  Freedom  (DOF).  This 
quantity  is  defined  as  number  of  variables  being  fit  subtracted  from  number  of  data  points 
used. 


New  estimates  of  each  parameter  were  tested  at  each  iteration.  The  new  estimates 
were  created  and  applied  using  the  bisection  method.  The  original  estimate  was  adjusted 
in  decreasing  intervals  described  by  Equation  3.14. 
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Estimate *ew  =  Estimate ±(0.5')  ■  Estimatex0rigional 


(3.14) 


where  X represents  the  variable  being  estimated  and  i  is  the  iteration  step.  The  values 
above  and  below  the  initial  estimate  were  evaluated  at  each  iteration.  The  RMSE  for  the 
original  estimate.  The  estimate  above  and  the  estimate  below  were  calculated.  Due  to 
the  form  of  the  data,  F  was  the  first  variable  fit  at  each  iteration.  The  estimated  value  of 
Tbomb  was  used  with  the  three  estimates  of  F.  The  F  estimate  resulting  in  the  minimum 
RMSE  was  then  taken  as  the  new  value  to  be  used  through  the  next  iteration.  The 
established  value  of  F  was  then  taken  with  the  three  estimates  of  bomb  temperature 
within  the  same  iteration  step.  The  value  of  Tbomb  resulting  in  the  lowest  RMSE  became 
the  value  used  in  the  next  iteration.  In  successive  iterations,  the  interval  of  the  variable 
change  is  halved  until  an  established  criteria  is  met.  The  criteria  used  for  this  research 
was  to  stop  iterating  when  the  value  of  the  step  size  fell  below  the  computer  resolution. 

When  fitting  two  variables,  the  process  of  adjusting  one  of  the  variables  affects 
the  best-fit  value  of  the  other  variable.  To  ensure  best  fit  values  were  as  accurate  as 
possible,  the  iteration  process  was  repeated.  The  new  iteration  set  was  started  using  the 
variable  values  established  in  the  previous  iteration  set.  This  process  was  repeated  until 
the  RMSE  for  the  final  estimated  values  remained  the  same  over  two  consecutive 
iteration  sets.  The  program  Table  Curve  was  used  to  verify  the  values  of  the  two 
variables.  The  values  found  by  TableCurve  were  within  1  -2  percent  of  the  values  found 
with  the  bisection  method  explained  above.  Details  of  the  TableCurve  comparison  can  be 
found  in  Chapters  IV  and  V. 
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IV.  Procedures 


Introduction 

Chapters  I,  II  and  III  discussed  the  theory  and  development  of  the  concepts  and 
tools  used  during  this  research.  This  chapter  describes  the  specific  details  of  the  process 
as  it  was  developed.  First,  an  explanation  of  the  format  of  the  ASCII  data  files  containing 
event  calibrated  scene  spectral  radiance  will  be  discussed.  This  will  lead  to  a  discussion 
on  how  to  determine  dimensions  of  the  ASCII  data  and  import  the  data  into  the 
FORTRAN  program.  Chapter  IV  continues  with  the  creation  of  PLEXUS  transmittance 
files  based  on  the  form  of  the  data  and  the  atmospheric  conditions  at  data  collection  time. 
Next  is  a  discussion  of  how  the  PLEXUS  transmittance  profile  is  used  to  remove 
atmospheric  signal  loss  from  the  data.  The  discussion  continues  with  an  explanation  of 
how  the  average  background  signal  was  determined  and  subtracted  from  the  data.  The 
process  used  to  fit  this  final  form  of  the  adjusted  data  and  the  statistical  quantities 
generated  is  discussed  next.  The  discussion  continues  with  a  description  of  the  data  files 
created  and  their  format.  Finally,  how  TableCutve  was  used  to  verify  curve  fit  results 
and  Mathematica  was  linked  to  the  FORTRAN  program  for  graphical  interpretation  of 
the  new  data  will  be  discussed. 

Data  Manipulation 

The  initial  approach  taken  in  this  effort  was  to  use  Mathematica  as  the  data 
manipulation  and  calculation  program.  Mathematica  provides  a  combination  of 
computational  and  graphical  tools.  However,  the  format  and  dimensions  of  the  data  files 
didn’t  work  efficiently  in  the  Mathematica  environment.  In  order  to  ingest  ASCII  data, 
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Mathematica  requires  the  dimensions  (number  of  rows  and  columns)  as  an  input. 
However,  the  data  files  had  different  dimensions  depending  on  sensor  used,  spectral 
resolution  selected  and  number  of  measurements  taken  during  a  data  collection  sequence. 
In  order  to  process  data  files  without  calculating  dimensions  first,  FORTRAN  was 
selected  as  the  program  to  be  used  for  data  manipulation  and  iterative  calculations.  The 
final  version  of  the  FORTRAN  code  is  included  as  Appendix  1 . 

The  process  of  creating  a  FORTRAN  program  began  with  the  following 
objectives.  First  the  program  should  be  able  to  read  the  data  regardless  of  dimensions. 
Second,  the  program  needed  to  be  able  to  place  the  data  into  arrays  to  be  used  during 
calculations.  Third,  the  code  had  to  be  able  to  combine  different  types  of  data  since  the 
input  data  did  not  always  have  the  same  dimensions.  Fourth,  the  program  had  to  accept 
and  use  different  parameters.  Finally,  the  program  had  to  create  new  data  files  containing 
the  results  of  the  calculations. 

The  inputs  to  the  program  would  be  the  locations  of  the  raw  data  created  in 
previous  research,  the  location  of  the  PLEXUS  transmittance  profile,  the  names  and 
locations  of  the  data  files  created  by  the  program  and  all  user  variables  used  in  the 
calculations.  The  final  version  of  the  FORTRAN  program  performed  all  functions  as 
directed. 

All  of  the  Input/Output  (10)  information  was  placed  in  a  file  named  “infile.txt”. 
The  10  file  could  be  changed  by  editing  the  file  directly  or  by  running  a  Mathematica 
script  built  to  edit  the  10  file.  The  Mathematica  code  used  to  create  the  10  file  and 
provide  graphical  data  plots  is  included  as  Appendix  2. 
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The  files  containing  the  detonation  event  data  were  identified  and  organized.  One 
ASCII  data  file  represents  the  total  calculated  scene  spectral  radiance  at  all  time  steps 
recorded.  Each  data  file  is  organized  in  rows  representing  the  intensity  value  at  a  given 
frequency  for  all  time  steps.  Columns  represent  the  value  at  a  given  time  step  for  all 
frequencies.  Table  4.1  shows  the  format  of  a  sample  ASCII  data  file. 


Table  4.1  ASCII  Data  format  of  Detonation  Event 


Frequency  (cmM) 

t„=i 

t„=2 

tn 

tn=nt-l 

tn=nt 

499.527 

Datau 

Datai>2 

Data^n 

Datai  nt_i 

Datalnt 

501.455 

Data2,i 

Data2,2 

Data2,n 

Data2,nt-i 

Data2,n, 

503.3839 

Data3>i 

Data3)2 

Data3,n 

Data35nt-i 

Data3,nt 

freqm 

Datamj 

Datam,2 

Data^n 

Data^nt-i 

Datalnt 

5998.1771 

Datanw_i  i 

Datanw_ii2 

Datanw_i2 

Datanw-i,nt-i 

Datanw-i,nt 

6000.1058 

D  at  anw.  i 

Datanw,2 

DatanW;3 

Datanwnt-i 

DatanW;nt 

where  nw  represents  the  number  of  frequency  divisions  and  nt  is  the  number  if  time  steps 
in  each  data  set.  The  first  column  represents  the  frequency  for  each  row  (in 
wavenumbers).  The  remaining  columns  held  data  values  at  different  time  steps.  Since 
there  is  a  column  for  frequency  values  the  total  number  of  columns  in  the  raw  data  is 
nt+1  and  the  dimensions  of  the  data  is  nw  by  nt+1. 

It  is  important  to  recognize  that  the  frequency  values  are  spaced  at  equal  intervals 
but  are  not  integer  values.  This  is  considered  when  applying  the  PLEXUS  created  data. 
PLEXUS  data  is  created  with  frequency  intervals  at  integer  values. 
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Once  the  program  read  the  data,  the  information  could  then  be  processed  in  a 
logical  order.  The  data  was  generally  processed  column  by  column.  This  order 
represents  an  interrogation  of  all  frequencies  simultaneously  as  time  progressed. 


PLEXUS  Application 

The  PLEXUS/MODTRAN  program  was  discussed  in  Chapters  I  and  II.  The 
program  creates  a  transmittance  profile  representing  atmospheric  absorption  or  scattering 
signal  loss.  Frequency  intervals  are  determined  by  the  user  of  the  software.  However, 
the  frequency  values  can  only  be  defined  in  integer  intervals.  In  order  for  a  row-by-row 
calculation  of  the  radiance  data  with  the  transmittance  data,  the  intervals  and  number  of 
rows  had  to  match.  Therefore,  a  new  transmittance  profile  was  created  for  each  data  set. 
The  new  profile  was  built  by  interpolating  the  value  of  the  transmittance  from  the  closest 
two  values  surrounding  each  raw  data  frequency  location.  The  interpolation  was  done 
using  linear  interpolation  based  on  divided  differences.  The  assumption  of  a  constant 
slope  is  made  over  the  interval  (xj,x2)  with  the  corresponding  values  yj  andy?.  A  value 
for  y,  can  be  found  at  a  given  value  x,-.  This  is  true  assuming  a  constant  slope  within  the 
interval.  With  a  constant  slope,  the  divided  difference  over  the  interval  (xj,Xi)  is  the  same 
as  the  divided  difference  over  the  interval  (xi,x2)  as  shown  in  Equation  4. 1 

U~Xl)  _  CV2-X1)  (A  1  > 


Solving  fory,  yields  Equation  4.2 


v  _  to-yi )<*<-*!)  ,  v 
Si  ~  {x2-xx)  h si 


(4.2) 


To  minimize  interpolation  error,  PLEXUS  transmittance  files  are  created  with 
frequency  resolution  of  the  integer  number  closest  to  the  frequency  resolution  of  the  data 
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being  corrected.  After  a  transmittance  profile  is  created  in  array  form,  the  original  data 
can  be  divided  by  the  transmittance  value  at  each  frequency.  If  the  transmittance  profile 
exactly  predicted  the  amount  of  signal  lost  to  the  atmosphere  at  each  frequency,  the  result 
would  be  an  exact  representation  of  the  signal  from  the  source.  However,  the  prediction 
is  not  exact  and  errors  in  the  process  appear  in  the  corrected  data.  These  errors  are 
greatest  where  values  of  transmittance  are  low.  Figures  4. 1  through  4.4  show  the  raw 
data  and  PLEXUS  corrected  data  for  both  sensors  at  a  time  step  immediately  after  bomb 
initiation. 
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Figure  4*1  HgCdTe  Raw  Data 
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Figure  4.2  HgCdTe  Corrected  Data 
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Figures  4.1  through  4.4  also  show  the  frequency  range  of  each  sensor.  The 
HgCdTe  sensor  produced  a  profile  with  more  noise  but  had  a  greater  frequency  range. 
The  InSb  sensor  produced  data  with  much  less  noise  but  the  frequency  range  did  not 
extend  into  the  Long  Wave  InfraRed  (LWIR)  region  of  the  spectrum.  These  sensor 
characteristics  are  important  to  remember  to  ensure  the  best  data  set  is  created  during  the 
PLEXUS  adjustment. 

Background  Signal  Subtraction 

The  final  correction  made  to  the  data  was  to  remove  the  average  background 
signal  prior  to  detonation  from  the  data  set.  The  explanation  of  the  background  signal 
subtraction  is  detailed  in  Chapter  III.  Without  subtracting  the  average  signal  profile  due 
to  background  radiation,  the  curve-fitting  process  would  fit  to  an  Equation  including  two 
Planck  functions  instead  of  one.  Additionally,  subtracting  the  average  background 
profile  eliminates  the  need  to  find  the  calibration  constant,  K,  to  be  used  in  the  two 
Planckian  version  shown  in  Equation  3.5. 

The  average  background  profile  was  built  by  summing  the  pre-detonation 
PLEXUS  corrected  data  profiles  over  a  number  of  time  steps.  The  number  of  time  steps 
used  was  then  divided  into  these  totals.  Although  the  first  two  time  steps  of  the  data 
appeared  to  be  representative  of  the  background  signal,  the  third  through  twenty- seventh 
time  steps  were  used  to  avoid  any  background  average  contamination  due  to  collection 
process  initiation.  Figures  4.5  and  4.6  show  typical  average  background  signal  profiles 
for  each  detector. 
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Scene  Spectral  Radiance  Background  Average  Scene  Spectral  Radiance  Background  Average 

(W/cmA2  sr  cmM)  (W/cmA2 sr  cmA-1 ) 


Figure  4.5  Average  HgCdTe  Background  Signal 


Figure  4.6  Average  InSb  Background  Signal 
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With  the  average  background  profile  subtracted  from  the  data,  the  only  data 
remaining  in  each  time  step  prior  to  event  initiation  is  noise.  No  information  could  be 
inferred  from  the  pre-detonation  data.  However,  the  conditions  filling  the  field  of  view 
are  known  quantities.  The  field  of  view  is  observing  the  background  which  is  basically  at 
the  ambient  atmospheric  temperature.  The  variable  F  representing  the  relative  fractional 
field  of  view  filled  by  the  event  is  zero  prior  to  detonation.  To  avoid  processing  the  data 
during  the  curve-fitting  process,  a  threshold  ‘trigger’  was  built  into  the  program  which 
determined  when  the  data  at  each  time  step  would  be  processed.  The  method  selected  was 
to  find  the  average  pre-background  subtracted  intensity  within  the  frequency  range  of 
2400-2800  wavenumbers.  This  interval  was  selected  for  three  reasons.  First,  the  interval 
is  included  in  the  range  of  both  sensors.  The  second  reason  this  region  was  selected  was 
the  PLEXUS  values  were  generally  greater  than  0.4.  The  final  reason  for  selecting  this 
interval  was  the  largest  increase  in  intensity  was  expected  in  this  region  for  the 
temperatures  of  the  events  expected.  The  average  background  signal  within  this  region 
was  easily  found  since  averages  of  the  entire  data  range  were  calculated  during 
background  average  calculation.  This  average  value  was  then  compared  to  the  average 
within  the  same  interval  of  the  actual  data  at  each  time  iteration. 

If  this  average  was  larger  than  the  actual  data,  the  time  step  was  given  values  of 
scene  temperature  equals  300  Kelvin  and  relative  fractional  field  of  view  due  to  event,  F, 
equals  zero.  The  first  time  the  average  value  fell  below  the  corrected  data  average,  event 
initiation  is  assumed  and  the  curve  fitting  process  begins  for  that  time  step.  The  values  of 
1500  Kelvin  for  temperature  and  2.0  for  F  were  used  as  initial  estimated  for  the  first  time 
step  processed.  In  successive  time  steps,  the  values  determined  during  the  fitting  process 
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of  the  prior  time  step  were  used  as  initial  estimates.  Once  both  corrections  were  applied 
to  the  data,  the  program  started  process  of  fitting  the  corrected  data  to  the  function. 

Curve  Fitting 

A  subroutine  was  built  within  the  program  to  test  the  estimates  being  made  during 
the  bisection  process  described  in  Chapter  III.  The  subroutine  created  statistical  values 
representing  how  well  the  function  matched  the  actual  data  using  the  estimates  provided. 
These  values  were  used  during  the  bisection  method  to  converge  to  the  best  values  of 
Tbomb  and  F.  The  form  of  these  statistical  values  are  shown  in  Equations  4.3  through  4.9. 

nw  ^ 

Sum  of  Squares  due  to  Error:  SSE  =  ^(S**  -  Sj,Tb°mh )  (4.3) 

M 

(Sum  of  Residuals  Squared) 


Sum  of  Squares  about  Mean: 

nw  f 

SSM  =  £  sj-s 

J= 1  v  J 

(4.4) 

Coefficient  of  Determination: 

(4.5) 

Degrees  of  Freedom 

DOF =n-m 

(4.6) 

Mean  Square  Error 

MSE=^r 

(4.7) 

Fit  Standard  Error  (Se) 

RMSE  =  yjMSE 

(4.8) 

F-statistic 

T?  -  MSR 
p  -  JTSF 

(4.9) 

Sj  represents  the  corrected  data  at  each  frequency.  Sj’Tbomb  is  the  value  at  each  frequency 
given  by  the  function  using  the  estimates  of  Tbomb  and  F .  S  is  the  mean  value  over  all 
frequencies  in  the  time  step.  The  total  number  of  data  points  used  is  n  and  m  is  the 
number  of  variables  in  the  function  being  fit.  The  exact  form  of  these  statistics  was  taken 
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from  the  TableCurve  3-D  Users  Manual  (23,E-5).  The  form  was  duplicated  since 
TableCurve  was  used  to  verify  fit  results  derived  by  the  program  and  show  graphical 
output. 

Data  Created 

An  execution  of  the  program  would  create  eight  files  in  ASCII  format.  The  first 
output  file  created  by  the  program  was  a  PLEXUS  transmittance  profile  with  interpolated 
frequency  values  matching  the  frequencies  of  the  event  data  file.  The  next  two  files  were 
created  by  solving  for  temperature  using  Planck’s  Equation.  These  files  were  created  to 
investigate  the  composition  of  the  main  data  file.  The  second  file  created  was  a 
temperature  as  a  function  of  frequency  and  raw  intensity  at  a  selected  time  step.  The 
third  data  file  built  was  a  temperature  file  for  all  times  and  frequencies.  The  fourth 
output  data  file  contained  all  the  statistical  fit  values  described  above  for  each  time  step 
processed.  The  fifth  ASCII  file  contained  the  profile  of  the  average  background  signal 
prior  to  detonation.  The  sixth  file  held  the  average  value  of  the  2400  through  2800  wave 
number  frequency  interval.  The  seventh  file  was  the  event  temperature  found  by  the 
curve  fitting  process  as  a  function  of  time.  The  eighth  output  file  contained  the  value  of 
the  variable  F  at  each  time  step  determined  by  the  fit.  The  final  data  file  created  during 
program  execution  was  the  final  value  of  RMSE  describing  the  fit  as  a  function  of  time. 

Interpretation  of  New  Data 

Reducing  the  original  data  was  accomplished  by  a  FORTRAN  program  created  as 
part  of  this  research.  The  FORTRAN  program  could  be  initiated  by  a  Mathematica  script 
which  provided  the  10  details.  Data  generated  through  execution  of  the  program  was 
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interpreted  using  two  additional  software  programs.  The  first  of  these  two  programs  was 
TableCurve  2-D  Version  5.0.  TableCurve  was  used  to  check  the  fit  of  the  derived 
variables  to  the  data  using  the  formula  derived  in  Chapter  III.  TableCurve  also  proved 
valuable  in  providing  a  graphical  representation  of  Planckian  curve  overlaying  the  data 
population.  The  second  program  used  to  interpret  the  derived  data  was  Microcal  Origin 
Version  6.0.  Microcal  Origin  organizes  data  in  spreadsheet  format.  The  spreadsheet 
format  proved  useful  in  displaying  three-dimensional  graphs  used  to  compare  data. 
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V.  Results 


Introduction 

This  chapter  details  the  results  of  the  data  processing  and  analysis  accomplished 
during  this  research  effort.  The  results  indicate  event  characterization  may  be  possible 
after  data  reduction.  Furthermore,  the  large  data  files  containing  the  event  could  be 
quickly  reduced  with  a  few  manipulations.  The  original  ASCII  data  files  were  between 
3-10  Mega-Bytes  each.  Each  file  was  reduced  with  the  program  described  in  previous 
chapters  into  three  manageable  20  Kilo-Byte  files  containing  the  temporal  description  of 
relative  fractional  field  of  view,  F,  bomb  temperature  and  fit  error.  Due  to  the  large 
amounts  of  data,  only  a  portion  of  the  events  and  data  will  be  described  in  this  chapter 
and  in  the  following  appendices.  For  continuity,  the  events  covered  in  this  work  are  the 
same  as  the  events  investigated  in  previous  research. 

Discussion  will  start  with  a  general  overview  of  the  three  dimensions  of  the  initial 
data.  This  is  followed  with  a  discussion  of  the  data  processing  as  it  was  done  step-by 
step.  This  discussion  includes  the  processes  used  to  verify  each  step.  Undesirable  effects 
on  the  data  from  the  processing  are  described.  Finally,  qualitative  examination  of  the 
graybody  characteristics  and  the  features  versus  time  curves  are  discussed. 

Data  Overview 

The  creation  of  data  manipulation  and  calculation  tools  was  a  success.  The 
programs  created  reduce  large  data  files  consisting  of  abstract  quantities  of  wavenumber 
and  scene  spectral  radiance  representing  the  IR  signal  of  an  event.  The  large  data  sets  are 
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converted  to  understandable  data  representing  the  temporal  behavior  of  the  event.  The 
initial  form  of  the  data  is  shown  in  Figure  5.1. 


Figure  5.1  Three-Dimensional  View  of  a  Typical  Event.  (17) 

This  event  was  collected  at  16  cm'1  spectral  resolution  using  the  InSb  detector. 
Each  wavenumber  bin  is  roughly  8  cm'1  apart.  Each  scan  represents  a  0.047  second  time 
step.  This  figure  represents  only  1.88  seconds  of  event  data.  The  entire  data  file  included 
27  seconds  of  data.  Vertical  axis  is  calibrated  scene  spectral  radiance  as  calculated  by  the 
MR- 154  and  Acquire  software.  Atmospheric  absorption  causes  the  zero  values  located 
between  intensity  peaks. 

A  single  slice  of  the  Figure  above  displays  the  spectral  profile  at  a  single  time  step 
measurement.  Figure  5.2  shows  the  spectral  profile  of  a  time  step  immediately  after 
detonation. 
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The  spectral  resolution  for  the  figure  is  4  cm1.  The  interval  of  the  frequency  steps 
for  both  detectors  is  1 .929  cm'1.  The  spectra  collected  with  the  InSb  detector  is 
represented  by  the  darker  line  ranging  from  1800  to  6000  cm'1.  The  HgCdTe  spectra  is 
the  more  noisy  data  ranging  from  500  to  6000  cm'1.  Atmospheric  absorption  bands 
degrade  or  eliminate  the  signal  in  several  frequency  intervals.  The  primary  atmospheric 
absorption  molecules  in  these  intervals  are  superimposed  on  the  FIGURE. 

The  noise  levels  in  the  HgCdTe  detector  increase  with  increasing  wavenumber. 
The  maximum  noise  range  for  the  HgCdTe  is  +/-  0.0001  W/(cm2  str  cm1).  The 
corresponding  noise  range  for  the  InSb  detector  is  approximately  +/-  0.00002  W/(cm2  str 
cm'1).  Spectra  collected  with  the  InSb  detector  contained  less  noise;  therefore,  the  data 
most  often  analyzed  was  collected  with  the  InSb  detector. 
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The  events  processed  include  data  collected  at  two  different  spectral  resolutions,  4 
cm1  and  16  cm1.  The  specific  spectral  features  of  the  events  are  not  investigated  in  this 
research.  The  work  done  here  eliminates  the  spectral  information  of  the  data  and  infers 
the  temporal  behavior  of  the  event’s  size  temperature  and  graybody  fit  error.  However, 
the  spectral  resolution  possibly  influenced  the  results  of  the  data  processing  sequence. 
This  possibility  is  discussed  later.  During  data  processing,  the  program  could  be 
instructed  to  include  instantaneous  profiles  of  the  data  at  a  given  time  step.  This 
capability  was  included  in  the  program  so  further  research  of  data  created  by  the  MR- 154 
could  be  investigated  through  its  temporal  or  spectral  characteristics.  The  specific  data 
manipulations  applied  to  the  data  can  be  explained  using  this  feature.  This  will  be  done 
by  following  the  data  shown  in  Figure  5.2  through  each  step  of  the  program  as  it 
executes. 

After  the  event  data  file  is  identified,  a  transmittance  profile  must  be  created  using 
the  atmospheric  conditions  during  data  collection  as  PLEXUS  input.  The  resolution  of 
the  created  transmittance  file  should  be  at  the  closest  integer  value  to  the  spectral 
resolution.  For  this  data  at  1.929  cm'1,  a  transmittance  file  is  created  at  2  cm'1  intervals. 
The  first  data  manipulation  is  to  interpolate  the  transmittance  values  to  values  at 
frequency  steps  in  the  data  file.  Once  this  is  accomplished,  the  transmittance  values  are 
divided  into  the  original  data.  Frequency  locations  with  transmittance  values  lower  than 
a  cutoff  value  provided  by  the  user  are  set  to  zero.  The  resultant  data  after  this 
manipulation  is  shown  in  Figure  5.3 
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Figure  5.3  PLEXUS  Adjusted  Scene  Spectral  Radiance 

Figure  5.3  shows  the  results  of  atmospheric  adjustment.  Division  by  cutoff  values 
too  large  causes  a  loss  of  data,  cutoff  values  too  low  result  in  overestimated  values.  The 
value  of  0.4  was  used  in  all  transmittance  corrections.  The  average  values  at  each 
frequency  prior  to  detonation  are  then  calculated  to  build  a  background  scene  spectral 
radiance  profile.  This  average  pre-detonation  signal  is  subtracted  as  described  in  Chapter 
IV  to  correct  the  calibration  near  ambient  temperatures.  The  average  background  profile 
subtracted  fr  om  the  post  detonation  data  were  shown  in  Figures  4.5  and  4.6.  Since  the 
magnitudes  of  the  average  background  profiles  were  orders  of  magnitude  smaller  than  the 
detonation  data,  the  final  form  of  the  data  is  effectively  unchanged. 

The  final  manipulation  of  the  adjusted  data  is  curve  fitting.  The  form  of  the 
function  to  describe  the  event  as  derived  in  Equation  3.6.  This  Equation  is  recalled  as 
Equation  5.1. 
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The  corrected  data  is  fit  to  this  Equation  as  described  in  Chapter  IV.  At  the  end  of 
each  time  iteration,  the  final  values  of  F,  Tbomb  and  the  RMSE  are  recorded.  These 
parameters  were  compared  to  the  values  computed  by  TableCurve.  Figures  5.4  and  5.5 
show  the  final  form  of  the  data  and  the  best-fit  curve  found  by  TableCurve  for  the  InSb 
and  HgCdTe  sensors  respectively. 


rA2=0. 45909276  DF  Adj  rA2=0. 45791 943  FitStdErr=0‘.00030773896  Fstat=783.3924 


a=0. 88483663 
b=  158 1.4343 


Figure  5.4  InSb  TableCurve  Data  Verification 
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rA2=0.57559026  DF  Adj  rA2=0. 57480938  FitStdErr=0.00040234865  Fstat=  1475.5604 
a=0. 98030393 
b=1577.5432 


Figure  5.5  HgCdTe  TableCurve  Data  Verification 


Notice  the  offset  of  the  curve  from  the  heavily  populated  data  sections.  This 
effect  was  considered  when  selecting  a  cutoff  for  the  transmittance  value.  When  a  lower 
value  was  used,  much  of  the  data  in  the  750  to  1250  wavenumber  region  was  lost.  This 
would  sometimes  cause  erroneous  curve  fitting  during  the  decay  phase  of  the  event 
evolution.  Chapter  VI  includes  a  discussion  on  possible  corrections  to  this  problem  as 
part  of  further  research. 

During  program  execution,  a  file  was  created  containing  the  final  values  of  the 
variables  and  the  statistical  quantities  describing  the  quality  of  the  fit  of  these  variables. 
As  a  comparison,  the  values  of  the  time  step  followed  above  are  compared  using  values 
from  this  file.  Table  5.1  shows  the  values  of  the  fit  variables  and  RMSE  determined  by 
the  program  compared  to  the  same  values  from  TableCurve. 
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Table  5.1  Variable  Fit  Comparison 


InSb 

F 

InSb 

Tbomb 

InSb 

RMSE 

HgCdTe 

F 

HgCdTe 

Tbomb 

HgCdTe 

RMSE 

FORTRAN 

0.89847 

1575.070 

3.0780e-4 

0.97714 

1578.887 

4.0235e-4 

TableCurve 

0.88484 

1581.434 

3.0774e-4 

0.98030 

1577.543 

4.0235e-4 

Overall,  the  values  calculated  by  the  FORTRAN  program  were  veiy  close  to  the 
values  derived  by  TableCurve.  A  statistical  analysis  of  an  entire  data  file  was  performed 
to  test  the  program  through  a  detonation  sequence.  The  data  files  used  included  the  time 
step  evolution  described  above. 

As  part  of  its  fitting  procedures,  TableCurve  calculates  values  above  and  below 
each  variable  representing  the  95%  confidence  interval  for  each  variable.  In  most  cases, 
the  FORTRAN  program  calculated  values  within  the  95%  confidence  intervals.  The  only 
time  the  values  strayed  from  the  TableCurve  intervals  was  with  the  InSb  during  the  decay 
phase  of  an  event.  This  is  believed  to  be  due  to  the  lack  of  data  below  1800 
wavenumbers.  Without  the  low  wavenumber  reference,  the  data  is  missing  a  reference 
for  half  of  the  Planckian  curve.  The  95%  confidence  intervals  calculated  are  included  in 
Figures  5.6  through  5.9  showing  the  InSb  and  HgCdTe  fit  comparison  of  variables  for  a 
Radiant  Brass  3B  Event  number  3. 
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-Relative  Fractional  Field  of  View 


Figure  5.8  Event  3  HgCdTe  Temperature  Comparison 


Figure  5.9  Event  3  HgCdTe  Fractional  Field  of  View  Comparison 


The  deviation  from  the  TableCurve  values  shown  in  Figure  5.7  is  probably  caused 
by  the  frequency  range  limitations  of  the  InSb  as  discussed  above.  The  deviation  of  our 
values  of  F  from  the  TableCurve  values  occurs  as  the  intensity  of  the  collected  signal  is 
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decreasing  rapidly  during  event  decay.  Figures  5.8  and  5.9  show  that  the  noise  on  the 
HgCdTe  detector  causes  the  confidence  interval  to  be  wider.  However,  the  FORTRAN 
code  generated  values  fit  the  TableCurve  values  for  all  data  from  the  HgCdTe  detector. 

Event  Classification 

During  his  work,  Orson  integrated  the  intensity  value  at  all  frequencies  for  each 
time  step  (17).  In  this  manner,  an  energy  versus  time  profile  is  created.  He  estimated  this 
profile  described  the  event  as  a  decaying  gray  body.  The  results  of  this  method  indicated 
two  different  time  evolution  patterns.  The  two  time  modes  were  identified  to  classify 
different  types  of  explosives.  Figures  5.10  and  5.11  from  his  thesis  show  these  two 
modes: 


Figure  5.10  Evolution  Mode  One  (17)  Figure  5.11  Evolution  Mode  Two  (17) 

Figure  5.4  shows  the  first  time  evolution  mode.  This  mode  was  the  expected  time 
decay  mode  and  has  two  well-defined  peaks.  This  pattern  described  the  phases  of 
initiation,  afterburn  and  decay.  The  second  time  evolution  mode  showed  an  initiation 
followed  by  a  growth  period  resulting  in  the  afterburn  peak  becoming  the  largest  apparent 
radiance  value.  The  decay  rate  of  this  mode  was  normally  slower  than  mode  one. 
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The  data  produced  from  this  work  shows  the  two  modes  also  can  be  used  to 
describe  the  temporal  evolution  of  the  event  Temperature.  Using  data  from  events  of  the 
same  explosion  species  produced  the  same  time  mode  profile  when  data  collection 
conditions  were  similar.  Some  collection  events  showed  inconclusive  results  due  to 
variability  in  look  angle,  impact  surface  and  impact  vector.  Data  adversely  affected  by 
these  factors  generally  displayed  the  same  temporal  profile  during  the  first  two  seconds 
following  detonation.  The  convention  of  identifying  events  as  explosive  types  A  and  B  is 
consistent  with  prior  research  in  this  area.  Two  species  in  particular  could  be  described 
by  different  time  modes. 

The  first  event  species  investigated  here  is  the  large  type  B  explosive.  Large  type 
B  detonations  were  clearly  identifiable  by  distinct  afterburn  feature  in  the  temperature 
profile  at  approximately  0.5  seconds  after  bomb  initiation.  Figures  5.12  through  5.14 
show  the  temporal  profiles  of  two  large  B  events,  collected  at  different  spectral 
resolutions. 


- Event#  15  HgCdTe 

16cm'1 

- Event#l  5  InSb 

16cm'1 

— I - . - 1 - . - 1 - 1 - 1 

2  3  4  5 

Time  (sec) 


Figure  5.12  Large  B  Temperature  Profiles 
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Figure  5.13  Large  B  Fractional  Field  of  View  Profiles 


Figure  5.14  Large  B  RMSE  Profiles 


Data  collected  on  large  B  events  consistently  displayed  the  afterburn  feature  from 
0.75  to  1.00  seconds  in  the  temperature  profile.  The  afterburn  is  followed  by  a  steady 
decay  out  to  5  seconds.  The  expected  slow  rise  of  the  variable  F  is  observed  clearly  in 
higher  frequency  resolution  data. 

The  second  type  of  data  investigated  was  collected  on  small  events  using  type  A 
explosives.  This  event  represents  the  largest  number  of  data  sets  in  ASCII  format. 
Figures  5.15  through  5.17  show  the  temporal  profiles  of  small  A  events. 
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Figure  5.15  Small  A  Temperature  Profiles 


Time  (sec) 


Time  (sec) 


Figure  5.16  Small  A  Fractional  FIELD  OF  VIEW  Profiles 
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small  A  detonations  displayed  an  instantaneous  initiation  followed  by  rapid  decay 
beginning  almost  immediately.  The  decay  rate  decreased  toward  the  end  of  the  event 
sequence.  It  may  be  possible  to  fit  these  temporal  profiles  to  an  exponential  decay 
model.  The  fractional  field  of  view  versus  time  profiles  were  highly  variable.  The 
general  pattern  seemed  to  be  a  smaller  and  nearly  constant  size  throughout  the  first  half 
of  the  event  sequence.  During  the  second  half  of  the  event  sequence,  about  2  seconds 
after  detonation,  the  size  of  the  event  indicates  a  gradual  growth  event  size. 

Data  from  other  species  was  analyzed  as  well.  Graphical  representations  of  these 
species  were  consistent  with  explosive  type.  Small  A  explosives  displayed  the  same  time 
evolution  features  as  large  A  and  medium  A  detonations.  Similarly,  large  B  explosives 
displayed  the  same  time  evolution  features  as  small  B  and  medium  B  detonations. 

The  FORTRAN  and  Mathematica  programs  created  should  be  able  to  read  and 
interpret  any  BOMEM  MR- 154  generated  ASCII  data  set  with  little  or  no  modification. 
Future  research  in  this  area  should  be  much  easier  with  the  programs  created  during  this 
thesis. 
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VI.  Conclusions 


Introduction 

This  chapter  addresses  the  conclusions  reached  through  the  process  of 
manipulating  data  collected  during  previous  research.  The  source  of  the  data  investigated 
was  the  calibrated  detonation  event  data  collected  using  the  Bomem  MR- 154  FTIR  to 
capture  detonation  signatures  during  the  Radiant  Brass  tests  at  Fallon  NAS.  The  chapter 
will  begin  with  a  short  discussion  on  the  process  used  and  the  data  created.  This  will  be 
followed  with  the  distinguishing  features  of  the  data  created.  Finally,  a  short  discussion 
of  the  overall  results  and  suggestions  for  future  research  are  presented. 

Summary 

A  FORTRAN  program  and  Mathematica  notebook  were  created  and  merged  to  a 
single  tool  that  was  used  to  manipulate  data  and  build  new  data  sets.  The  final  version  of 
these  tools  proved  to  be  efficient  in  handling  large  data  files  of  different  sizes  and 
processing  the  data  they  contained.  New  data  files  were  created  for  each  data  set 
processed.  The  new  data  files  included  temporal  profiles  of  event  temperature,  relative 
fractional  FOV,  RMSE  for  data  processed  at  each  time  step,  statistical  values  for  each 
time  step  processed  and  intensity  files  representing  profiles  of  an  individual  time  step 
selected  by  the  user.  The  code  was  created  so  future  research  efforts  could  manipulate 
the  data  differently  with  minimal  changes  to  the  code.  The  program  was  used  to  process 
over  thirty  data  files  in  the  course  of  this  research. 
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Conclusions 


The  investigation  of  the  derived  data  show  identifiable  features  for  different 
explosive  species.  Some  data  did  not  clearly  display  the  anticipated  features  after 
processing.  These  anomalies  appeared  to  be  from  differences  in  the  parameters 
surrounding  the  event  such  as  observation  elevation,  angle  of  attack  and  impact  surface. 

In  most  of  these  cases  the  features  were  as  expected  but  were  not  definite. 

The  resulting  data  created  verified  the  description  of  detonation  events  as 
decaying  graybodies.  The  data  from  large  type  B  explosives  consistently  produced  a 
clear  afterburn  feature.  The  fractional  field  of  view  for  all  type  B  events  indicated  a  peak 
value  corresponding  with  event  initiation  followed  immediately  by  a  distinct  drop-off  and 
steady  value.  Type  A  explosives  showed  an  exponential  decay  of  temperature  from  the 
peak  value  of  initiation.  The  fractional  field  of  view  display  for  type  A  explosives 
showed  a  steady  profile  without  an  initial  peak.  The  form  of  the  RMSE  generally 
showed  the  same  time  mode  as  its  associated  temperature  profile.  These  observations 
indicate  that  event  species  can  possibly  be  determined  by  temperature  and/or  fractional 
field  of  view  profiles. 

Discussion 

The  findings  of  this  research  indicate  that  the  concept  of  reducing  large 
spectroradiometic  data  files  into  data  representing  event  temporal  features  may  be  a 
viable  means  of  event  classification.  Future  research  could  refine  the  transmittance 
profile  creation  and  the  curve-fitting  procedures  used  to  better  accomplish  these 
reductions. 
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The  temporal  evolution  of  the  fractional  field  of  view  profiles  may  be  verifiable 
using  IR  imagery  corresponding  to  each  event.  IR  imagery  was  collected  during  the 
Radiant  Brass  data  acquisition.  Comparing  this  IR  imagery  with  the  fractional  field  of 
view  profiles  may  verify  the  results  of  this  research  effort. 

Collection  of  a  wider  range  of  explosive  events  would  provide  a  better  data  set. 
Completion  of  analysis  on  all  data  already  collected  would  also  provide  a  better  analysis. 
The  low  angle  of  observation  seemed  to  be  the  major  factor  in  degradation  of  initial  data 
confidence.  An  elevated  viewing  angle  during  future  data  collection  opportunities  may 
increase  the  value  of  data  collected  during  future  research. 

The  process  of  reversing  the  atmospheric  absorption  could  possibly  be 
accomplished  much  more  accurately.  The  transmittance  values  created  with  PLEXUS 
appeared  to  overestimate  the  absorption  of  signal  at  certain  frequencies.  The  product  of 
this  over  estimation  was  transmissivity  values  lower  than  actual  conditions.  This  caused 
overestimated  values  of  adjusted  scene  spectral  radiance  from  the  atmospheric  correction 
process.  The  PLEXUS  software  used  user  input  of  the  parameters  ambient  temperature, 
path  length,  geographic  location,  time/date  of  measurement  and  path  elevation  as  inputs 
to  the  MODTRAN  model.  Relative  humidity,  atmospheric  pressure  and  percentage  of 
CO2  are  inferred  from  other  inputs.  A  better  initiation  of  MODTRAN  using  actual  values 
of  pressure,  humidity  and  C02  content  would  almost  certainly  provide  a  better 
atmospheric  correction  profile. 

Finally,  the  curve-fitting  process  used  in  the  FORTRAN  program  created  could  be 
improved.  During  program  execution,  the  curve  fitting  process  takes  1-4  seconds  for 
each  time  step.  Over  500  iterations  were  required  to  reach  the  final  variable  values. 
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TableCurve  found  better  results  in  less  than  20  iterations  with  a  process  taking  tenths  of  a 
second  on  average.  Better  numerical  method  techniques  would  almost  certainly  provide 
better  results. 
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Appendix  A:  FORTRAN  Analysis  Code 


Capt  William  F.  Bagby 

input  and  interoperation  of  ASCII  data  files 

6  Nov  00 

calculating ,  ver  4.05 


MODULE  global_data 

IMPLICIT  NONE 
SAVE 

REAL,  PARAMETER ::  h=6.626E-34 
REAL,  PARAMETER ::  c=2.9979E10 
REAL,  PARAMETER ::  k=1.38E-23 
REAL,  PARAMETER ::  C1=1.2E-12 
REAL,  PARAMETER ::  C2=1.439 
END  MODULE  global_data 
MODULE  sqerror 

USE  global_data 
IMPLICIT  NONE 
CONTAINS 

SUBROUTINE 


!  Planck's  constant! 

!  Speed  of  Light  in  cm/sec 
!  Boatman’s  constant 
!  2*h*cA2 
!  h*c/k 


LEASTERROR(s  1  ,bkav,i  ,fgss,tgss,nw,sqdev,rsq,fstat,mse,msr  ,ssm,rmean,sse,n  1 ) 


IMPLICIT  NONE 
INTEGER,  INTENT(IN) ::  i 
REAL,  INTENT(IN) ::  sl(:,:) 
REAL,  INTENT (IN) ::  bkav(:,:) 
REAL,  INTENT(IN) ::  fgss 
REAL,  INTENT(IN) ::  rmean 
REAL,  INTENT (IN) ::  tgss 
REAL,  INTENT(OUT) ::  sqdev 
REAL,  INTENT(OUT) ::  rsq 
REAL,  INTENT (OUT) ::  Fstat 
REAL,  INTENT (OUT) ::  MSE 
REAL,  INTENT (OUT) ::  MSR 
REAL,  INTENT(OUT) ::  SSM 
REAL,  INTENT (OUT) ::  SSE 
REAL ::  DOF 
REAL ::  RMSE 
REAL ::  ssll,ss22,ss33 
REAL  ::  wave, wave  1 
REAL  ::  diff,diffsum,diffsq 
REAL  ::  meandiff,meandiffsq 
REAL  ::  diffsumsq 
INTEGER  ::n,nl 
INTEGER ::  m 
INTEGER ::  j,nw 


!  Time  step  being  evaluated 
!  Intensities  array 
!  avg  background  intensity 
!  Estimated  F  value 
!  Mean  intensity  vale  at  time  step  i 
!  Estimated  Temperature  value 
!  Root  mean  square  error 
!  rA2  (coefficient  of  determination) 

!  F  statistic  (MSR/MSE) 

!  Mean  Square  Error 
!  Mean  Square  Regression 
!  Sum  of  Squares  about  Mean 
!  Sum  of  Squares  due  to  Error 
!  Degrees  of  Freedom  n-m 
!  Root  Mean  Square  Error  (Fit  Std  Error) 

!  Planckian  temporary  values 
!  Input  values  of  frequency  and  radiance 
!  Temporary  stats  values 
!  Temporary  stats  values 
!  Temporary  stats  values 
!  Number  of  data  points  used 
!  Number  of  variables  to  be  fit  (2) 

!  counters  used  to  step  through  frequencies 


!  PRINT*, "i^  ",i,"  tguess=  ",tgss,"  fguess=  ",fgss  ! ! !  ,"T  bomb=  M,tsave  !  record  fguess  at 

time  step  i 

!  PRINT*, "nw=",nw,M  rmean-',  rmean 

m=2  !  Number  of  variables  to  be  fit  (2) 

n=0 

diffsum=0 

diffsq=0 

diffsumsq=0 
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sqdev=0 

MSE=0 

RMSE=0 

MSR=0 

Fstat=0 

meandiff=0 

meandiffsq=0 

SSE=0 

SSM=0 

DO  j=l,nw 

wave=s  10,1)  !  Wave  number 

!  wavel=PLEXUS  corrected  intensity-background  intensity 

wavel=slO,i) 

IF  (wavel  .ne.  0)  THEN 

wave  1 =s  1 0  ,i)-bkav0 ,2) 

ssll=(fgss)  *  (Cl  *  wave**3/(Exp(C2  *  wave/tgss)-l) ) 
ss22=ssll  !  plankian  radiance  value  for  T=tempguess 
ss33=wavel  !  value  from  PLEXUS  corrected  data 
n=n+l 

meandiff=  ss33  -  rmean  !  difference  from  mean  squared 
ELSE  IF  (wavel  ==  0)  THEN 
ss22=0 
ss33=0 

meandiff=  0  !  difference  from  mean  squared 

END  IF 

diff=  ss22  -  ss33  !  difference  between  radiance  values 

difFsq=diff*diff  !  to  be  used  as  absolute  RMS  error 

SSE=SSE+diffsq  !  cumulative  SSE  errors. 

meandiffsq=meandiff*meandifF  !  diff  from  mean  squared  for 

absolute  RMS  error 

SSM=SSM+meandiffsq  !  cumulative  SSE  errors. 

END  DO 

DOF=n-m 

nl=n 

MSE=SSE/DOF 

RMSE=SQRT(MSE) 

sqdev=RMSE 

MSR=(SSM-SSE)/(m-l)  !  Mean  Square  Regression 

Fstat=MSR/MSE  !  F-statistic 

rsq=l-SSE/SSM 

END  SUBROUTINE  LEASTERROR 
END  MODULE  sqerror 

Program  inputdata 
USE  global_data 
USE  sqerror 

IMPLICIT  NONE 

INTEGER.PARAMETER ::  in3=3,inl=l,in2=2 
INTEGER.PARAMETER ::  out=10,outl=ll,out2=12 
INTEGER.PARAMETER ::  out3=13,out4=14,out5=15 
INTEGER, PARAMETER ::  out6=16,out7=17,out8=18 
INTEGER.PARAMETER ::  out9=  19, out  10=20, out  11=21 
INTEGER.PARAMETER ::  outl2=22,outl3=23,outl4=19 
INTEGER ::  iwn.it 
counter 


!  unit  numbers  for  reading  external  files 
!  unit  numbers  for  writing  external  files 
!  unit  numbers  for  writing  external  files 
!  unit  numbers  for  writing  external  files 
!  unit  numbers  for  writing  external  files 
!  unit  numbers  for  writing  external  files 
!  wavenumber  and  time  step  iteration 
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INTEGER ::  nw,nwl,nt,nn 
INTEGER ::  ios,iosl,ioread 
INTEGER ::  commas, crs 
INTEGER ::  n,ij,jl,jj jjl 
INTEGER ::  m,ml,ntl,z,FF,ii 
INTEGER ::  error 
INTEGER ::  length, flag2,flagl 
INTEGER  ::  add  1, flag, tflag,Fflag 
INTEGER ::  iter,iterl,nnn 

REAL,ALLOCATABLE  ::  s(:,:) 

REAL, ALLOCAT ABLE  ::  sl(:,:) 

REAL, ALLOCAT ABLE  ::  T(:,:) 
REAL,ALLOCATABLE  ::  Tl(:,:) 

REAL, ALLOCAT  ABLE  ::  TT(:,:) 
REAL,ALLOCATABLE ::  FT(:,:) 
time  data 

REAL, ALLOCAT  ABLE  ::  times(:) 

REAL, ALLOCAT  ABLE  ::  wn(:) 

REAL, ALLOCAT  ABLE  ::  tm(:,:) 

REAL,  ALLOCAT  ABLE  ::  trnl(:,:) 

REAL, ALLOCAT  ABLE  ::  SSM1(:,:) 
REAL, ALLOCAT  ABLE  ::  bkav(:,:) 
REAL,ALLOCATABLE ::  ss2(:) 

REAL, ALLOCAT  ABLE  ::  ss3(:) 

REAL  ::  delt,delwn,deltp 
REAL ::  tll,tl2,sll,sl2,sl3,s4,s5 
REAL  ::  final wn 

REAL  ::  temps, tmpF, lower, lowers 
REAL ::  p,pl,avg,mm 
REAL  ::  tguess,fguess, bin I,bin2,bin3, test 
REAL  ::  abinl,abin2,abin3, wave, stempo, rad 


!  total  bins  (intensities,  PLEXUS,  time) 

!  datafile  read/write  error  flags 
!  total  number  of  commas  and  returns 
!  integer  counters  for  loops 
!  integer  counters  for  loops 
!  error  code  for  input 

!  length  of  character  string  from  first  record 
!  integers  for  calculating  bin  averages 
!  integers  for  calculating  bin  averages 

!  intensity  data  @  each  time  and  wave# 

!  PLEXUS  corrected  intensity  data 
!  Raw  data  temperature  array 
!  PLEXUS  data  temperature  array 
!  Temperature  vs.  time  data 
!  Event  fraction  of  FIELD  OF  VIEW  vs. 

!  time  steps  array 
!  wave  number  array 
!  transmittance  values  from  PLEXUS 
!  transmittance  interpolated  from  PLEXUS 
!  Array  of  errors  for  each  time  step  fit 
!  avg  background  subtracted 
!  holds  timestep  mean  value 
!  holds  timestep  bin  mean  value 
!  time  step  size,  wave  number  step  size 
!  temporary  value  holders  for  building  arrays 
!  calculated  final  wave  number  value 


REAL  ::  sse,time  !diffsum,diffsumtot,diffsumsq,diff,diffsq,  various  curve  fitting  reals 
REAL  ::  abinlold,abin2old,abin3old,oldtguess,fguesstart,fguessl ,mse,msr,ssm,rl,r lsum,rmean 
REAL  ::  diffsqold,sqdev,delF,oldsq,fgup,sqdn,sqdv,squp,fgdn,fstat,fitsq,fitsqold  !,fitsqdiff 
REAL ::  tguesstart,oldsql,tgdn,tgup,tguessl,oldsqt,sqdnt,squpt,rsq,rsqold,rsqtemp  !,rsqdiff 
REAL  ::  rsqup,rsqdn,oldrsq, ambient, tsave,sqtold,bkval,valbk,binav,abin, launch, FicF, flow 


CHARACTER ::  firstchar 
CHARACTER*  100000  ::  charstring 
CHARACTER  ::  char  string  1 
CHARACTER*  10  ::  stringvar 
CHARACTER  (LEN=1)::  tempchar 


!  temporary  char  to  determine  number  of  records 
!  temporary  char  to  determine  elements  per  record 
!  temporary  char  to  determine  datafile  size 
!  temporary  char  to  determine  datafile  size 
!  temporary  char  to  determine  datafile  size 


CHARACTER  (LEN=40) ::  infile,name,namel  !  input  data  files  path  and  name 
CHARACTER  (LEN=40) ::  outfile,outfilel,outfile2  !  names  output  files  created 

CHARACTER  (LEN=40) ::  outfile3,outfile4,outfile5  !  names  output  files  created 
CHARACTER  (LEN=40) ::  outfile6,outfile7,outfile8  !  names  output  files  created 
CHARACTER  (LEN=40) ::  outfile9,outfilel0,outfilel  1  !  names  output  files  created 

CHARACTER  (LEN=40) ::  outfilel2,outfilel3,outfilel4  !  names  output  files  created 


! ! ! ! !  PROMPT  for  and  receive  name  and  path  of  ASCII  data  file  to  be  used 
infile- d:\infile.txt' 

OPEN  (UNIT=2,  FILE=infile,STATUS="OLDM,IOSTAT=ios) 

READ  (UNIT=2,FMT='(A40),,IOSTAT=ios),  name  !  Location  of  ASCII  Data  File  IN 

READ  (UNIT=2,FMT-(A40)',IOSTAT=ios),  namel  !  Location  of  ASCII  Transmittance  file 

READ  (UNIT=2,FMT-(A40),,IOSTAT==ios),  outfile2  !  Interpolated  Transmittance  file 


READ  (UNIT=2,FMT='(A40)',IOSTAT=ios),  outfilel 

READ  (UNIT=2,FMT='(A40)',IOSTAT=ios),  outfile 

READ  (UNIT=2,FMT=’(A40)',IOSTAT=ios),  outfile4 

READ  (UNIT=2,FMT='(A40)',IOSTAT=ios),  outfile5 

READ  (UNIT=2,FMT='(A40)',IOSTAT=ios),  outfile6 

READ  (UNIT=2,FMT='(A40)',IOSTAT=ios),  outfile7 

READ  (UNIT=2,FMT='(A40)',IOSTAT=ios),  outfile8 

READ  (UNIT=2,FMT='(A40)',IOSTAT=ios),  outfile9 

READ  (UNIT=2,FMT=*,IOSTAT=ios),  delt 

READ  (UNIT=2,FMT=*,IOSTAT=ios),  jj 

READ  (UNIT =2 ,FMT  =  * ,10 ST AT=ios) ,  lower 

READ  (UNIT=2,FMT='(A40)',IOSTAT=ios),  outfflelO 

READ  (UNIT=2,FMT=’(A40)',IOSTAT=ios),  outfilel  1 

READ  (UNIT=2,FMT='(A40)',IOSTAT=ios),  outfilel2 

name=TRIM(name) 

name  1  =TRIM(name  1 ) 

outfile2=TRIM(outfile2) 

outfilel=TRIM(outfilel) 

outfile3=TRIM(outfile3) 

outfile=TRIM(outfile) 

outfile4=TRIM(outfile4) 

outfile5=TRIM(outfile5) 

outfile6=TRIM(outfile6) 

outfile7=TRIM(outfile7) 

outfile8=TRIM(outfile8) 

outfile9=TRIM(outfile9) 

outfile  1 0=TRIM(outfil  e  1 0) 

outfilel  l=TRIM(outfilel  1) 

outfilel  2=TRIM(outfilel  2) 


!  Intensity  at  selected  timestep 
!  Interpolated  PLEXUS  applied  data 
!  Temps  at  selected  timestep 
!  Temps  for  all  timesteps 
!  Avg  Temp  vs.  time 
!  F  vs.  Time 
!  SSM  vs.  time 
!  Binvalue  vs.  time 
!  Timestep  size 
!  Selected  timestep  location 
!  cutoff  value  for  PLEXUS 
!  Average  Background  profile  (outfilel  0) 
!  adjusted  data  1  time  step 
!  fit  errors  data  listing; 


CLOSE  (UNIT=2) 

PRINT  *,  "ASCII  INPUT  FILE;  (name)=  ",name 

PRINT  *,  "PLEXUS  INPUT  FILE;  (namel)=  ",namel 

PRINT  *,  "Interpolated  PLEXUS  File;  (outfile2)=  ",outfile2 

PRINT  *,  "Intensities  vs.  wavenumber  for  1  timestep  file;  (outfilel)=  ",outfilel 

PRINT  *,  "Output  ASCII  divided  by  PLEXUS  tansmittance  file;  (outfile)=  ", outfile 

PRINT  *,  "Output  temperature  vs.  wave  number  file  for  1  timestep;  (outfile4)=  ",outfile4 

PRINT  *,  "Output  ASCII  Temperature  Data  file;  (outfile5)=  ",outfile5 

PRINT  *,  "Output  Temperature  vs.  time  file;  (outfile6)=  ",outfile6 

PRINT  *,  "Output  F  vs.  time  file;  (outfile7)=  ",outfile7 

PRINT  *,  "fit  error  vs.  time  file;  (outfile8)=  ",outfile8 

PRINT  *,  "2400-2800  cmA-l  avg  amplitude  vs.  time  file;  (outfile9)=  ",outfile9 

PRINT  *,  "Average  Background  profile  (outfilel 0)=  ",outfilel0 

PRINT  *,  "PLEXUS  and  Background  adjusted  data  1  time  step  (outfilel  1)=  ".outfilel  1 

PRINT  *,  "fit  errors  data  listing;  (outfilel2)=  ",outfilel2 

PRINT  *,  "Temporal  Resolution,  the  size  of  the  timestep;  (delt)=  ”,delt 

PRINT  *,  "Location  for  selected  time  step  to  be  viewed=  ”,jj 

PRINT  *,  "Lower  PLEXUS  BOUNDS,  or  the  cutoff  for  which  wavenumbers  values  are  used=  ".lower 


!!  ******  OPEN  ASCII  INPUT  FILE 

OPEN  (UNIT=3,  FILE=name,STATUS="OLD",IOSTAT=ios) 
IF  (ios==0)  THEN 
GO  TO  44 

ENDIF 
PRINT  *, " " 


75 


PRINT  *  "***************************************************************  n 
PRINT  *,  "Unable  to  open  data  file,  checF  path  and  filename  and  try  again" 

PRINT  *,  "The  file  below  has  an  invalid  filename  or  path  " 

PRINT  *,  name 
PRINT  *, "  " 

PRINT  *,  "Press  return  to  continue.  " 

PRINT  *, " " 

READ* 

GO  TO  900 

44  CONTINUE 

crs=0 

commas=0 

n=0 

nw=0 

nt=0 

DO  n=  1,1 8888  !  find  dimensions  of  array 

crs=crs+l 

IF  (crs  /=  1)  THEN  !  counts  the  number  of  records  in  the  data  set 

READ  (UNIT=3,FMT='(Al)',IOSTAT=ios),  firstchar 
SELECT  CASE  (ios)  !  Test  for  end  of  file  and  exit  loop 

CASE  (:-l)  !  end  of  file,  no  more  data 

GO  TO  100 

CASE  (1 :)  !  error  reading  ASCII  data 

error  =  -2 

PRINT  *,  "ERROR  Reading  ASCII  DATA!!!" 

GO  TO  100 

END  SELECT 

ELSE  IF  (crs  ==  1)  THEN  !  Test  for  commas  seperating  data  entries  and  count 

READ  (UNIT=3  ,FMT=’(A 1 00000)’),  charstring  !  Read  1st  record 

length=LEN(TRIM(charstring))  !  eliminate  excess  blanF  spaces  in  charstring 
DO  i=l , length  !  calculate  length  of  trimmed  char  string 

IF  (charstring(i:i)  ==  ',')  THEN  !  Count  commas 

commas=commas+l 

END  IF 
END  DO 

END  IF 
END  DO 
100  Continue 
nt=commas+l 
nw=crs-l 
PRINT  *, "  " 

PRINT  *,  "ASCII  file  wavenumber  bins,  nw=”,nw  !  number  of  records  or  wavenumber  bins 

PRINT  *,  ”  ” 

PRINT  *,  "ASCII  file  time  step  columns,  nt=",nt  !  number  of  columns  or  time  steps 

ntl=nt-l  !  Use  nt-1;  first  entry  in  each  record  is  a  wave  number  not  intensity 

ALLOCATE  (s(nw,nt))  !  Allocate  size  of  intensities  array 

ALLOCATE  (times(ntl))  !  Allocate  size  of  time  step  array 

ALLOCATE  (wn(nw))  !  Allocate  size  of  wavenumber  array 

!!!!!!!!!!!!  Create  the  1-Dtimestep  array  !!!!!!!!! 

DO  i=l,ntl 

times(i)=(i- 1 )  *delt 
END  DO 

CLOSE  (UNIT=3) 


OPEN  (UNIT=3,  FILE=name,STATUS="OLD",POSITION="REWIND",IOSTAT=ios) 
!!!!!!!!!!!  Create  the  2-D  intensities  array  !!!!!!! 
wn=0 
DO  i=l,nw 

READ  (UNIT=3  ,FMT=*,IOSTAT=ioread)(s(i  j)  j= 1  ,nt) 
wn(i)=s(i,l) 

IF  (ios/=0)  THEN 

PRINT  *,  '("Error  ", 15, "while  reading  ",A)',  name 

ENDIF 
END  DO 

del wn=wn(2)-wn(  1 ) 
final  wn=delwn*nw 

!!!!  OPEN  PLEXUS  INPUT  FILE 

OPEN  (UNIT=1,  FILE=namel,STATUS="OLD",IOSTAT=ios) 

IF  (ios==0)  THEN 
GO  TO  45 

ENDIF 
PRINT  *, "  " 

PRINT  *  ”***************************************************************  " 
PRINT  *,  "Unable  to  open  data  file,  checF  path  and  filename  and  try  again" 

PRINT  *,  "The  file  below  has  an  invalid  filename  or  path  " 

PRINT  *,  namel 

PRINT  *,  "Press  return  to  continue.  " 

READ* 

GO  TO  900 
45  CONTINUE 

crs=0 

n=0 

nwl=0 

DO  n=  1,9888  !  Outer  loop  counts  the  number  of  records  in  the  PLEXUS  data  set 

crs=crs+l 
IF  (crs  /=  0)  THEN 

READ  (UNIT=  1  ,FMT='(A1  )',IOSTAT=ios),  firstchar 
SELECT  CASE  (ios)  !  Test  for  end  of  file  and  exit  loop 

CASE  (:-l)  !  end  of  file,  no  more  data 

GO  TO  110 

CASE  (1:)  !  error  reading 

error  =  -2 

PRINT  *,  "ERROR  Reading  PLEXUS  DATA!!!" 

GO  TO  110 

END  SELECT 

ENDIF 
END  DO 
110  CONTINUE 
nwl=crs-l 
PRINT  V  " 

PRINT  *,  "PLEXUS  wavenumber  bins,  nwl=",nwl 

ALLOCATE  (tm(nwl,2))  !  Allocate  size  of  origional  PLEXUS  array 

ALLOCATE  (trn  1  (nw,2))  !  Allocate  size  of  interpolated  PLEXUS  array 

CLOSE  (UNIT=1) 

!!!!!!!!!!!  Create  the  2-D  PLEXUS  array  !!!!!!! 

OPEN  (UNIT=  1 ,  FILE=namel  ,STATUS="OLD",POSITION="REWIND",IOSTAT=ios) 
DO  i=l,nwl 
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READ  (UNIT=  1  ,FMT=*,IOSTAT=ioread)(trn(ij)j=l  ,2) 

IF  (ios/=0)  THEN 

PRINT  *,  "Error  while  reading ",  namel 
GO  TO  111 

ENDIF 
END  DO 
111  CONTINUE 
PRINT  *, "  " 

CLOSE  (UNIT=1) 

! ! !  Create  a  new  transmittance  file  with  same  wavenumber  bins  as  the  ASCII  file  (interpolated) 

OPEN  (UNIT=2,FILE=outfile2, ST ATUS= "REPLACE", IOSTAT=iosl)  !  Output  interpolated 

PLEXUS  outfile2 

OPEN  (UNIT=  1 1  ,FILE=outfilel  ,STATUS="REPLACE",IOSTAT=iosl)  !  Selected  step  ASCII 

Data  outfilel 

trnl=0 

i=l 

124  CONTINUE 
IF  (i  >  nw)  THEN 
GOTO  123 

ENDIF 

sll=wn(i) 

sl2=s(ijj) 

WRITE  (UNIT=21,FMT=*)  si  1  ,sl2  !  interpolated  PLEXUS  file 

DO  m=l,nwl 

IF  (wn(i)>=  tm(m,l))  THEN 
ml=m+l 

IF  (wn(i)  <  tm(ml,l))  THEN 
tll=wn(i) 

1 1 2=(wn(i)-trn(m,  1 ))  *  ((trn(m  1 ,2)-tm(m,2))/  (tm(m  1 , 1  )-tm(m,  1  )))+trn(m,2) 

tml(i,l)=tll 

tml(i,2)=tl2 

WRITE  (UNIT=21,FMT=*)  trn  1  (i,  1  ),trn  1  (i,2)  !  Create  interpolated  File 
i=i+l 

GO  TO  124 

ELSE  IF  (wn(i)  >=  tm(ml,l))  THEN 
CONTINUE 

ENDIF 

ELSE  IF  (wn(i)  <  tm(m,l))  THEN 
PRINT*,  "I  was  here" 

GO  TO  123 

ENDIF 
END  DO 
123  CONTINUE 
ENDFILE  21 
CLOSE  (UNIT=21) 

ENDFILE  2 
CLOSE  (UNIT=2) 

! ! ! !  Divide  the  intensity  by  the  transmittance  to  approximate  original  emittance 
PLEXUS  adjusted  outfile 
ALLOCATE  (sl(nw,nt)) 

DO  i=l,nw 

sl(i,l)=wn(i) 

DO  j=2,nt 

IF  (tml(i,2)  <  lower)  THEN 
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tl  1=0 
sl(ij)=tll 

ELSE  IF  (tml(i,2)  >=  lower)  THEN 
tl2=l  /  tml(i,2) 
tll=tl2*s(ij) 
sl(ij)=tll 

END  IF 
END  DO 
END  DO 

ALLOCATE  (Tl  (nw,nt))  !  array  of  all  temps  derived  from  Plexus  corrected  data 

ALLOCATE  (T(nw,2))  !  array  of  temps  at  selected  time  step 

OPEN  (UNIT=4,FILE=outfile4,STATUS="REPLACE",IOSTAT=iosl)  !outfile4=Derived  temps  for  1 

time  step 

lowers=.0000001 

DO  i=l,nw  ! !  Create  a  derived  temperature  array 
tll=wn(i) 
tl(i,l)=tll 
DO  FF=l,nt 

IF  (ABS(sl(i,FF))  <  lowers)  THEN 
tl  1=0 

tl(i,FF)=tll 

ELSE  IF  (ABS(sl(i,FF))  >=  lowers)  THEN 
IF  (sl(i,FF)  >=  0.0)  THEN 
s4=ABS(sl(i,FF)) 
tl  l=(c  *  h  *  sl(i,l))  /  (k) 
tl2=  log(  ((2*  h  *  c**2  *  (sl(i,l))**3)  /(s4)  )+l  ) 
flagl=1.0 

ELSE  IF  (sl(i,FF)  <  0.0)  THEN 
s4=ABS(sl(i,FF)) 
tl  l=(c  *  h  *  sl(i,l))  /  (k) 
tl2=  log(  ((2*  h  *  c**2  *  (sl(i,l))**3)  /(s4)  )+l ) 
flagl=-1.0 

END  IF 

tl(i,FF)=(flagl  *(tl  l/tl2)) 

END  IF 
END  DO 

WRITE  (UNIT=4,FMT=*)  wn(i),tl(ijj)  !  Temperatures  from  PLEXUS  Corrected  data 
END  DO 
ENDFILE  4 
CLOSE  (UNIT=4) 

OPEN  (UNIT=5,FILE=outfile5,STATUS="REPLACE",IOSTAT=iosl)!  Out  ASCII  Temperature 

Data=outfile5 

DO  i=l,nw 

WRITE  (UNIT=5,FMT=*)  (tl(i,j)  j=l,nt)  !  ASCII  Temperature  vs.  wave#  (all  times) 

ENDDO 
ENDFILE  5 
CLOSE  (UNIT=5) 

CLOSE  (UNIT=3) 

!!!!!  CURVE  FIT  THE  DATA  TO  MATCH  A  TEMPERATURE  AND  COEFFICIENT  !!!!!!! 
ALLOCATE  (bkav(nw,2)) 

ALLOCATE  (ss2(ntl-l)) 

ALLOCATE  (ss3(ntl-l)) 

ALLOCATE  (TT(ntl-l,2)) 

ALLOCATE  (FT(ntl-l,2)) 

ALLOCATE  (SSMl(ntl-l,2)) 
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bkval=0 

valbk=0 

FicF=.0000001 

binav=0 

tguess=1500 

fguess=2 

tflag=0 

abinlold=10 

ambient=300 

flow=.02 

iter=0 

fguessl=fguess 

tguessl=tguess 

n=0 

nn=0 

bkav=0 

abin=0 

!  outfile  1 2="c  :\errors.txt" 

OPEN  (UNIT=12,FILE=outfilel2,STATUS="REPLACE",IOSTAT=iosl)!  Out  ASCII  Temperature 
Data=outfile5 

WRITE  (UNIT=  1 2,FMT=*)  "lower  PLEXUS  bounds=  ".lower 

WRITE  (UNIT  =  1 2,FMT=  * )  "ASCII  INPUT  FILE;  ".name 

WRITE  (UNIT=12,FMT=*)  "PLEXUS  INPUT  FILE;  ".namel 
WRITE  (UNIT=12,FMT=*)  "Interpolated  PLEXUS  File;  ",outfile2 
WRITE  (UNIT=12,FMT=*)  "ASCII  divided  by  PLEXUS  file;  ".outfile 
WRITE  (UNIT=12,FMT=*)  "ASCII  Temperature  Data  file;  ",outfile5 
WRITE  (UNIT=  1 2,FMT=*)  "Temperature  vs.  time  file;  ",outfile6 
WRITE  (UNIT=12,FMT=*)  "Fractional  FIELD  OF  VIEW  vs.  time  file;  ",outfile7 

WRITE  (UNIT=12,FMT=*)  "fit  error  vs.  time  file;  ",outfile8 
WRITE  (UNIT=12,FMT=*)  "2400-2800  cmA-l  avg  vs.  time;  ",outfile9 
WRITE  (UNIT=  1 2,FMT=*)  "Temporal  Resolution,  ",delt 

WRITE  (UNIT=12,FMT=*) "  " 

DO  i=2,ntl 

rmean=0 
rlsum=0 
add 1=0 
abinl=0 
bin  1=0 
!  n=0 

DO  j=l,nw 

rl=sl(j,2) 
rlsum=rlsum+  rl 

IF  (slO'.i)  -ne.  0)  THEN  !  Bin  1  2400-2800  cmM 

addl=addl+l 
IF  (slOM)  >  2400)  THEN 

IF  (sl(j,l)<  2800)  THEN 
binl=binl  +  sl(j,i) 
n=n+l 

END  IF  !  (si (j.l)  <  2800) 

END  IF  !  (sl(j,l)  >  2400)  THEN 

END  IF 

END  DO  !  j=l,nw 
rmean=rlsum  /  addl 
abinl=binl  /  addl 
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ss2(i-l)=abinl 


!  average  signal  in  2400-2800  range  at  timestep  i 
!  average  signal  at  timestep  i 


ss3(i-l)=rmean 
IF  (i  >=  3)  THEN 

IF  (i  <  12)  THEN 

binav=abin  1  +binav 
nn=nn+l 

END  IF 

END  IF 

END  DO!  i=2,ntl 
binav=binav/nn 
bkav=0 
nnn=0 
DO  j=l,nw 

nnn=0 
DO  i=3,12 

nnn=nnn+1 

IF  (sl(j,i)  -ne.  0)  THEN 

bkav(j  ,2)=s  1  (j  ,i)+bkav(j  ,2) 

ELSE  IF  (sl(j,i)  ==  0)  THEN 
bkav(j,2)=0 

END  IF 
END  DO 
END  DO 

Print*, "addl=",addl 
Print*,"nnn= "  ,nnn 
DO  j=l,nw 

bkav(j,l)=sl(j,l) 

IF  (bkav0',2)  .ne.  0)  THEN 

bkav(j  ,2)=bkav(j  ,2)/nnn 

END  IF 
END  DO 
launch=0 

OPEN  (UNIT=1  l,FILE=outfilel  l,STATUS="REPLACE",IOSTAT=iosl) 

OPEN  (UNIT=  1 0,FILE=outfile,STATUS="REPLACE",IOSTAT=iosl) 

DO  i=l,nw 

bkav(i,2)=bkav(i,2)  !/addl 

bin l=sl(ijj)  -  bkav(i,2) 

WRITE  (UNIT=10,FMT=*)  sl(i,l),binl  !  PLEXUS  And  Bkgd  Corrected  1  time  step  jj 
WRITE  (UNIT=  1 1  ,FMT=*)  (bkav(i  j)  j=  1 ,2)  !  av  bkd  signal  profile 

ENDDO 
ENDFILE  11 
CLOSE  (UNIT=11) 

ENDFILE  10 
CLOSE  (UNIT=10) 

DO  i=2,ntl 

Fflag=0 

iter=0 

rsqold=l 

sqtold=0 

n=0 

rmean=ss3(i-l) 

test=ss2(i-l) 

IF  (test  <=  2.0  *  binav)  THEN  !  check  for  a  large  increase  in  intensity 
tguessl=ambient 
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fguess  1=0.0 
launch=0 
GO  TO  446 

ELSE  IF  (test  >  2.0  *  binav)  THEN !  check  for  a  large  increase  in  intensity 
IF  (launch  ==0)  THEN 
tguessl=1500 
fguess  1=2 
launch=launch  +  1 

END  IF 

442  Continue 

fguesstart=fguess  1 
tguesstart=tguessl 
iter=0 

Fflag=Fflag+l 

444  Continue 

flag=0 

443  Continue 

iter=iter+l 

IF  (iter  ==  120)  THEN 

PRint*,  "I  got  out  at  line  527" 

GO  TO  446 

END  IF 

delF=(.5  **  iter)  *  fguesstart  Ifind  MSE  at  fguess 

CALL 

LEASTERROR(sl  ,bkav,i,fguessl  ,tguessl  ,nw,sqdev,rsq,fstat,mse,msr,ssm,rmean,sse,n) 
oldsq=sqdev 
oldrsq=rsq 

fgdn=fguess  1  -delF  Ifind  MSE  below  fguess 

CALL 

LEASTERROR(sl,bkav,i,fgdn,tguessl,nw,sqdev,rsq,fstat,mse,msr,ssm,rmean,sse,n) 

sqdn=sqdev 

rsqdn=rsq 

fgup=fguess  1  +delF  Ifind  MSE  above  fguess 

CALL 

LEASTERROR(sl,bkav,i,fgup,tguessl,nw,sqdev,rsq,fstat,mse,msr,ssm,rmean,sse,n) 

squp=sqdev 

rsqup=rsq 

IF  (oldsq  <=  squp)  THEN 

IF  (oldsq  <=  sqdn)  THEN 
fitsqold=oldsq 
GOTO  451 

ELSE  IF  (oldsq  >  sqdn)  THEN 
IF  (squp  >  sqdn)  THEN 
fguess  l=fgdn 
fitsqold=sqdn 
GOTO  451 

END  IF 

END  IF 

ELSE  IF  (oldsq  >  squp)  THEN 
fguessl=fgup 
fitsqold=squp 
GOTO  451 

END  IF 
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451  CONTINUE 

deltp=(.5  **  iter)  *  tguesstart 
CALL 

LEASTERROR(sl,bkav,i,fguessl,tguessl,nw,sqdev,rsq,fstat,mse,msr,ssm,rmean,sse,n) 

oldsqt=sqdev 

rsqtemp=rsq 

tgdn=tguess  1  -deltp  Ifind  MSE  above  Tguess 

CALL 

LEASTERROR(sl,bkav,i,fguessl,tgdn,nw,sqdev,rsq,fstat,mse,msr,ssm,rmean,sse,n) 

sqdnt=sqdev 

rsqdn=rsq 

tgup=tguess  1  +deltp  Ifind  MSE  below  Tguess 

CALL 

LEASTERROR(sl  , bka  v,i,fguessl,tgup,nw, sqdev, rsq,fstat,mse,msr,ssm,rmcan,sse,n) 
squpt=sqdev 
rsqup=rsq 

IF  (oldsqT  <=  squpT)  THEN 

IF  (oldsqT  <=  sqdnT)  THEN  !  sqdev  is  less  than  up  or  down  step 

IF  (rsqold  ==  rsqtemp)  THEN 
IF  (flag  >=  90)  THEN 
GO  TO  446 

ELSE  IF  (flag  <  90)  THEN 
Flag=flag+1 
rsqold=rsqtemp 
go  to  443 

END  IF 

ELSE  IF  (rsqold  .ne.  rsqtemp)  THEN 
rsqold=rsqtemp 
IF  (1  <=  Fflag)  THEN 

IF  (Fflag  <=  75)  THEN 
GO  TO  442 

END  IF 

END  IF 
GO  TO  444 

END  IF 

ELSE  IF  (oldsqT  >  sqdnT)  THEN  !  sqdev  is  less  than  up  or  down  step 

IF  (squpt  >  sqdnt)  THEN 
tguess  l=tgdn 
rsqold=rsqdn 
IF  (1  <=  Fflag)  THEN 

IF  (Fflag  <=  75)  THEN 
GO  TO  442 

END  IF 

END  IF 
GO  TO  444 

END  IF 

END  IF 

ELSE  IF  (oldsqt  >  squpt)  THEN 
tguess  l=t  gup 
rsqold=rsqup 
IF  (1  <=  Fflag)  THEN 
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IF  (Fflag  <=  75)  THEN 
GO  TO  442 

END  IF 

END  IF 
GO  TO  444 

END  IF 

END  IF  !  (ss2(i-l)  <=  1.5  *  binav) 

446  CONTINUE 

fguess=fguessl 

tguess=tguessl 

fitsq=oldsqt 

time=i*delt 

If  (launch  .ne.  0)  THEN 
Print*,"" 

CALL 

LEASTERROR(sl,bkav,i,fguess,tguess,nw,sqdev,rsq,fstat,mse,msr,ssm,rmean,sse,n) 

Print*, "step  i=",i,"  at  time=",time,"  after", iter,"  iterations" 

PRINT*,"fguess=  ",fguess,"  tguess=  ",tguess  !!!  ,"T  bomb=  ",tsave  !  record  fguess  at 

time  step  i 

PRINT*,"SSE=",SSE,"  SSM=",SSM,"  nw=",nw 
PRINT*,"MSR=",MSR,"  MSE=",MSE,"  RMSE=",sqdev 
PRINT*, "Fstat=",fstat,"  rsq="qsq,"  n="qi 
WRITE  (UNIT=  1 2,FMT=*)"" 

WRITE  (UNIT=12,FMT=*)"step  i=",i,"  at  time=",time,"  after", iter,"  iterations" 
WRITE  (UNIT=12,FMT=*)"fguess=  ", fguess,"  tguess=  ",tguess  !!!  ,"Tbomb=  ",tsave 
!  record  fguess  at  time  step  i 

WRITE  (UNIT=12,FMT=*)”SSE=",SSE,"  SSM=",SSM,"  nw="qiw 
WRITE  (UNIT=12,FMT=*)"MSR=",MSR,"  MSE=",MSE,"  RMSE=",sqdev 
WRITE  (UNIT=12,FMT=*)"Fstat=",fstat,”  rsq=",rsq,"  n="qi 

END  IF 

TT(i-l,l)=times(i-l)  !  Create  timestep 

TT(i-l,2)=tguess  !  Temperature  converted  in  Kelvin 

FT(i-l,l)=times(i-l) 

FT(i-l,2)=fguess  !  FIELD  OF  VIEW  Fraction  guess 
SSMl(i-l,l)=times(i-l) 

SSMl(i-l,2)=sqdev  !  error  guess 

END  DO  !  i 

! ! !  Create  temperature  for  each  time  step  file 

OPEN  (UNIT=6,FILE=outfile6,STATUS="REPLACE",IOSTAT=iosl)  !  Out  ASCII  Temp  vs.  time 
=outfile6 

OPEN  (UNIT=7,FILE=outfile7,STATUS="REPLACE",IOSTAT=iosl)  !  Out  ASCII  F  vs.  time  =outfile7 
OPEN  (UNIT=8,FILE=outfile8,STATUS="REPLACE”,IOSTAT=iosl)  !  Out  ASCII  fit  errors 

DO  jj=  1  ,nt  1-1 

WRITE  (UNIT=6,FMT=*)  (TT(jjj)j=l,2)  !  Temperature  vs.  time  in  K 

WRITE  (UNIT=7,FMT=*)  (FT(jjj)j=  1 ,2)  !  FIELD  OF  VIEW  Fraction  guess  vs.  time 

WRITE  (UNIT=8,FMT=*)  (SSMl(jjj)j=l,2)  !  error  guess  vs.  time 

ENDDO  !jj 
ENDFILE  6 
CLOSE  (UNIT=6) 

ENDFILE  7 
CLOSE  (UNIT=7) 

ENDFILE  8 
CLOSE  (UNIT=8) 

ENDFILE  12 
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CLOSE  (UNIT=12) 

900  CONTINUE 
END  Program  inputdata 


Appendix  B:  Mathematica  Analysis  Code 

This  is  a  sample  Mathematica  notebook  used  to  create  an  10  file  for  the  FORTRAN 
code  and  then  execute  the  program.  Below  the  Run[“c:\prog.exe”]  command  are 
commands  developed  to  display  the  contents  of  the  data  files  built  during  program 
execution.  During  data  processing  two  of  these  notebook  scripts  were  created,  one  for 
the  A  sensor  on  the  instrument  and  the  other  for  the  B  sensor. 


«  Graphics' MultipleListPlot' 

h  =  6.626*  10A(-34);  c  =  2.9979*10^10;  k  =  1.38*10^-23); 
T1  =  300;  T2  =  873;  T3  =  1073;  T4  =  1173;  T5  =  1273; 


S[A_,  T_]  :=  ((10A16/[A]A5)*((2*h*cA2)/((Exp[(1000*c*h)/(  A*k*T)]  -  1))) 

(*  Frequency  Planck  curve  *) 

Sl[cr_,  T_]  :=  (2*h*cA2*[<x]A3)/((Exp[(c*h*[o-])/(k*T)]  - 1)) 

(*  Wave  Number  Planck  curve  *) 

TT[S_J<x_J  :=  (c*h*[cr])/(k*Log[(2*h*cA2*[cr]A3)/S  +  1])) 

(*  Inverse  Planck  to  get  Temperature  *) ))) 

lists  =  Table[ (j,  SID,  1100]},  {j,  2,  5990,  2}];  (*  1100K  Curve  *) 

(*  Transmissivity  file  *) 
tms  =  ReadList[var2a,  {Number,  Number}]; 

ListPlot[list5,  Axes  ->  True,  ImageSize  ->  400,  Frame  ->  True, 

FrameLabel  ->  {"Wavenumber  (cmA-l)",  "Radiance",  "For  1100  K",  "  "  }] 

ListPlotftrns,  ImageSize  ->  400,  Axes  ->  True,  Frame  ->  True, 

FrameLabel  ->  {"Wavenumber  (cmA-l)",  ’Transmittance",  "PLEXUS",  var2a  }] 

(*  300  K  curve  times  transmittance  from  plexus  *) 
compl  =  Tablet  {2*j  +  396,  (list5[[j  +  199,  2]]*trns[[j,  2]])},  {j,  2, 2800, 2}]; 

(*  Now  set  limits  for  same  wavenumber  outputs  *) 

ListPlot  [compl,  Axes  ->  True,  ImageSize  ->  400,  Frame  ->  True,  Plot  Joined  ->  True,  PlotRange  ->  { {-200, 
6000},  {-.00002,  .00078}},  FrameLabel  ->  {"Wavenumber  (cmA-l)",  "Spectral  Radiance  (W/cmA2  sr  cmA- 
1)", "  1100K",  ""}] 

(*  Define  Location  of  input  intensities  ASCII  File  to  be  used  *) 
varla  =  "d:Bigla.txt"; 

(*  Define  Location  of  input  PLEXUS  Transmittance  File  to  be  used  *) 
var2a=  "c:bigl.trn"; 

(*  Define  Location  for  output  interpolated  PLEXUS  Transmittance  File  *) 
var3a=  "d:TempNewtrn_a.tm"; 

(*  Define  Location  for  output  Temperatures  at  all  timesteps  data  *) 
var7a  =  "d:TempTcube_a.txt"; 
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(*  Define  Location  for  output  average  temperature  vs.  time  data  *) 
var8a=  "d:TempTemptime_a.txt"; 

(*  Define  Location  for  output  average  F  vs.  time  data  *) 
var9a  =  "d:TempFtime_a.txt"; 

(*  Define  Location  for  output  error  vs.  time  data  *) 
varlOa  =  "d:Temperrortime_a.txt"; 

(*  Define  Location  for  output  average  2400  -  2800  cmM  value  vs.  time  data  *) 
varl la  =  "d:Tempbinvaltime_a.txt"; 

(*  Define  time  step  size  in  seconds  *) 
delt  =  .0454; 

(*  Define  Location  for  fit  errors  and  statistics  *) 
varl4a  =  "diTempFitErrora.txt"; 

(*  Define  Location  for  output  Avg  bkgd  profile  *) 
varl2a  =  "d:Tempavbacka.txt"; 

(*  Define  data  step  number  to  be  used  for  instantaneous  values  *) 
step  =  41; 

(*  Define  Location  for  raw  output  intensity  at  1  time  step  *) 
var4a  =  "c:dataE98_03aRaw42a.txt"; 

(*  output  iterpolated  PLEXUS  applied  data  at  1  time  step  *) 
var5a  =  "c:dataE98_03aPlex42a.txt"; 

(*  PLEXUS  and  background  adjusted  data  at  1  timestep  *) 
varl3a=  "d:adjdata42a.txt"; 

(*  Define  Location  for  output  Temperatures  at  1  timestep  *) 
var6a  =  "d:Temps42a.txt"; 

(*  Define  lower  bounds  of  PLEXUS  to  be  used;  (Cutoff)  *) 
bounds  =  .4; 

varl  =  OutputForm[varla];  var2  =  OutputForm[var2a];  var3  =  OutputForm[var3a]; 
var4=OutputForm[var4a];  var5  =  OutputForm[var5a];  var6  =  OutputForm[var6a]; 
var7=OutputForm[var7a];  var8  =  OutputForm[var8a];  var9=OutputForm[var9a]; 
varlO=OutputForm[ varlOa];  varl  l=OutputForm[varl  la];  varl2  =  OutputForm[varl2a]; 
varl3=OutputForm[varl3a];  varl  4  =  OutputForm[varl4a]; 

Put[varl,  var2,  var3,  var4,  var5,  var6,  var7,  var8,  var9,  varlO,  varl  1,  delt,  step,  bounds,  varl2,  varl3, 
varl4,  "d:infile.txt"] 

Run["d:prog.exe"]; 
tstep  =  step*delt; 

(*  raw  data  at  one  timestep  *) 

data2  =  ReadList[var4a,  {Number,  Number}]; 

ListPlot[data2,  Axes  ->  True,  ImageSize  ->  400,  Frame  ->  True, 

PlotRange  ->  {{0,  6000},  {-.0001,  .0012}}, 

FrameLabel  ->  {"Wavenumber  (cmA-l)", 

"Scene  Spectral  Radiance  (W/cmA2  sr  cmA-l)",  "HgCdTe  Raw  Data",  ""}] 


(*  PLEXUS  Corrected  raw  data  at  one  time  step  *) 
corr  =  ReadList[var5a,  {Number,  Number}]; 

ListPlot[corr,  Axes  ->  True,  ImageSize  ->  400,  Frame  ->  True,  PlotRange  ->  { {0,  6000},  {-0.0001, 
.007}},  FrameLabel  ->  {"Wavenumber  (cmA-l)",  "Scene  Spectral  Radiance  (W/cmA2  sr  cmA-l)", 
"HgCdTe  Corrected  Data",  "  "  }] 

(*  Temperatures  from  raw  data  in  FORTRAN  *) 
temps  =  ReadList[var6a,  {Number,  Number}]; 

ListPlot[temps,  Axes  ->  True,  ImageSize  ->  400,  Frame  ->  True,  PlotRange  ->  {{0,  6000},  {-3000, 
3000}},  FrameLabel  ->  {"Wavenumber  (cmA-l)",  "Temperature  (Kelvin)",  var6a,  tstep}] 
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(*  Temperature  vs.  time  *) 

temps  =  ReadList[var8a,  {Number,  Number}]; 

ListPlot[temps,  Axes  ->  True,  Plot  Joined  ->  True,  ImageSize  ->  400,  Frame  ->  True,  PlotRange  ->  { {1, 
4},  {-100,  1650}},  FrameLabel  ->  {"Time  (sec)", "  Temperature  (K)",  var8a,  bounds  }] 

(*  K  vs.  time  *) 

temps  =  ReadList[var9a,  {Number,  Number}];  ListPlot[temps,  Axes  ->  True,  PlotJoined  ->  True, 
ImageSize  ->  400,  Frame  ->  True,  PlotRange  ->  {{1,4},  {-.1, 15}},  FrameLabel  ->  {"Time  (sec)",  "F", 
var9a,  bounds  }] 

(*  average  background  radiance  *) 

bkav  =  ReadList["c:bkav.txt",  {Number,  Number}]; 

ListPlot[bkav,  Axes  ->  True,  ImageSize  ->  400,  Frame  ->  True,  PlotRange  ->  {{000,  6050},  {-.0001, 
.0003}},  FrameLabel  ->  {"Time  (sec)",  "Scene  Spectral  Radiance  (W/cmA2  sr  cmA-l)",  "HgCdTe 
Average  Background  Signal",  ""  }] 

(*  2400  -  2800  cmA-l  avg  vs.  time  *) 
temps  =  ReadList[varlla,  {Number,  Number}]; 

ListPlot[temps,  Axes  ->  True,  PlotJoined  ->  True,  ImageSize  ->  400,  Frame  ->  True,  PlotRange  ->  { {- 
0.10,  11},  {-.001,  .004}},  FrameLabel  ->  {"2400-2800  cmA-l  avg",  "bin  average",  var  11a,  bounds  }] 


(*  error  vs.  time  *) 

temps  =  ReadList[var  10a,  {Number,  Number}]; 

ListPlot[temps,  Axes  ->  True,  PlotJoined  ->  True,  ImageSize  ->  400,  Frame  ->  True,  PlotRange  ->  { {- 
0.10,  11},  {-0.0000001,  .0009}},  FrameLabel  ->  {"2400-2800  cmA-l  avg",  "error",  varlOa,  bounds  }] 
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